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ABSTRACT 



Photometric rotational modulations due to starspots remain the most common 
and accessible way to study stellar activity. In the Kepler-era, there now exists precise, 
continuous photometry of ~150,000 stars presenting an unprecedented opportunity for 
Q_[ statistical analyses of these modulations. Modelling rotational modulations allows one 

to invert the observations into several basic parameters, such as the rotation period, 
spot coverage, stellar inclination and differe ntial rotat i on ra te. T he most widely used 
analytic model for this inversion comes from Budding (|1977l ) and iDorre n (1987), who 
considered circular, grey starspots for a linearly limb darkened star. In this work, 
we extend the model to be more suitable in the analysis of high precision photome- 
try, such as that by Kepler. Our new freely available Fortran code, macula, provides 
several improvements, such as non-linear limb darkening of the star and spot, a single- 
domain analytic function, partial derivatives for all input parameters, temporal partial 
derivatives, diluted light compensation, instrumental offset normalisations, differential 
rotation, starspot evolution and predictions of transit depth variations due to unoc- 
cultcd spots. Through numerical testing, we find that the inclusion of non-linear limb 
darkening means macula has a maximum photometric error an order-of-magnitudc 
less than that of iDorrenl ((T987), for Sun- like stars observed in the ifepZer-bandpass. 
The code executes three orders-of-magnitude faster than comparable numerical codes 
making it well-suited for inference problems. 



Key words: methods: analytical — techniques: photometric — stars: spots — plan- 
etary systems 



1 INTRODUCTION 
1.1 Stellar Activity 

A variety of cool stars with external convection envelopes 
have been observed t o exhibit magnetic activity similar t o 
that of the Sun (e.g. lKronlll947l : iMullanl Il97i IVogtlll975h . 
These magnetic fields, generated by cyclonic turbulence in 
the outer convection zone, penetrate the stellar atmospher e 
forming starspots, plages, networks, etc |Berdvug ina 2005). 
The study of these manifestations on other stars allows 
for crucial tests of stellar dynamo theory. For example, 
ISkumanic hi (|1972|) first suggested that rotation plays a key 
role in generating stellar activity. 

Since the discovery of rotationally modulated bright- 
ness variations due to starspots, photometry remains the 
most common technique for studying stellar activity. In 



particular, space-based photometric instruments have pro- 
vided man y high-cadence, precise light curves (e.g. HIP- 
PARCOS; Ivan Leeuwen et al.l Il997l ). Recently, the detec- 
tion of transiting extrasolar planets |Charbonneau et alj 
l2000l : iHenrv et al.l [2000) has led to a surge in the design 
and construction of precise photometric instruments. No- 
tably, the Kepler Mission has detected 2165 eclipsing bi- 
naries dSlawson et alj|201ll ) and 2321 planetary candidates 
l|Batalha et al.ll2012l ) with nearly continuous photometry at 
a precision of ~ 50ppm (for V ~ 12) per long-cadence ex- 
posure (29.4244 minutes). Preliminary analysis of the Kepler 
target stars (over 150,000) has revealed that those which ex- 
hibit periodic modula tion generally hav e a much higher am- 
plitude of variability (|Basri et alj|20lj ). One of the legacy 
products of the Kepler Mission will be a vast database of 
precise continuous photometry and the effective exploitation 
of this database will surely lead to deep insights into stellar 
activity. 
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1.2 Starspots 

Starspots are a very common source of photometric variabil- 
ity and have a diverse value to astronomers, varying from 
friend to foe. The presence of dark starspots on the surface of 
a rotating star induces periodic photometric variability due 
to the stellar rotation. An analysis of these rotational modu- 
lations allows for a determination of several basic properties 
of the star. The most accessible of these properties is the 
rotation period, which can often be inferred using a simple 
Lomb-Scargle periodogram or autocovariance analysis, and 
has several astrophysical uses. For example, the rotation pe- 
riod may b e used with gyrochronology to estimate the age 
of the star (Barnes 2009). Another e xamp le is demonstrated 
in the recent work of iHirano et al.l |2012h who show how a 
spectroscopic V sin/*, an estimate of the stellar radius (R,) 
and the rotation period allows one to infer the stellar incli- 
nation angle, I t . 

Employing spot-modelling codes allows for more infor- 
mation t han just the rotatio n period to be derived. For 
example, Wal ker et al.l ((2007) used rotational modulations 
alone to infer the stellar inclination angle for k 1 Ceti. Here 
the authors also showed how their measurement could be 
used to predict V sin J, and verified their solution was con- 
sistent with a spectroscopic determination. Further more, 
the authors were also able to estimate the differential rota- 
tion rate of k 1 Ceti, which was found to be reasonably close 
to Solar. 

If a star hosts a transiting planet which passes over 
a dark starspot, the transit light curve appear s to increase 
over the duration of the spot crossing event (e.g. lRabus et all 
2009). Detecting the same spot-crossing event in two con- 
secutive transits, which will have migrated along in longi- 
tude, allows the observer to infer a nearly coplanar spin- 
orbit angle. This technique has so far been successfully ap- 
plied to several ca ses including CoRoT - 2b ([Nutzman et aU 
l201ll) WASP-4b JSanchis-Oieda et all l201ll ) and Kepler- 
17b (|Desert et all 120111 ). In the case of WASP-4b, the re- 
sult was verified using the more traditional spec troscopic 
technique known as the Rossiter-McLaughlin effect (|Rossiterl 
1 1924 iMcLaughlinll 1924! ) . Modelling spot -crossing events re- 
mains outside of the scope of this work, but exploiting such 
phenomena is greatly aided by including information en- 
co ded in the out-o f -trans it photometry too, as pointed out 
bv lNutzman et all l|201ll) . 

In contrast to the examples given so far, other authors 
consider starspots to be a nuisance rather than a tool, due 
to their differing goals. For example, the "Hunt for Exo - 
moons with Kepler" (HEK) project (|Kipping et all 120121 ) 
anticipates that starspot crossings will be a source of false- 
positives for exomoon identification due to the morpho- 
logical similarities with planet-moon mutual events. Cross- 
referencing the transit anomaly with rotational modulations 
may be used to test whether the event is consistent with a 
starspot or not ( Sa nchis-Oieda et al.ll2012l ). 

Finally, even non-occulted starspots are a source of 
frustration in some arenas. In particular, these spots sub- 
tly change the perceived transit depth. Since spots vary in 
both time and wavelength, they are therefore capable of pro- 
ducing transit depth variations. This point is highly salient 
for those studying the atmospheres of exoplanets, who seek 
small chromatic variations in the depth. The study of ro- 



tational modulations can be used to correct th e res ulting 
transmission spectra, as shown in iDesert et all (|2009l ) for 
example. 

It is therefore clear that the study and interpretation of 
rotation modulations due to starspots is crucial to several 
areas of modern astronomical research. In this work, we aim 
to provide a revised model for such modulations which can 
account for several previously ignored effects. 

1.3 Numerical vs Analytic Models 

Inverting the brightness modulation of a star into a physi- 
cal map of the starspot coverage can be broached in several 
ways. The most successful technique is Doppler imaging aug- 
mented by precise photometry (e.g. see Vogt fc Penrod 1983; 
ICollier-Cameron et al.l Il994l ; iTuominen et al.l |2002| ). How- 
ever, using rotational modulation alone, as is the case for 
the Kepler Mission, is a more challenging problem. 

Eclipse mapping of eclipsing binaries (EBs) allows for 
greatly improved inversions of the spot coverage through 
photometry alone. This is usually done by numerically pixe- 
lating the star and applying t he maximum-entropy-metho d 
(MEM) to invert the map fe.g. lCollier-Cameron et al.lll997l ). 
In principle, a transiting planet offers the same opportunity 
as was discussed earlier for measuring spin-orbit alignments. 
However, eclipse-mapping using planets has several draw- 
backs. For example, the much smaller ratio-of-radii means 
only a thin-strip of the star is actually sampled. This means 
that inferences about the non-eclipsed portion of the star, 
which affect the perceived transit depth for example, must 
be made using the rotational modulations. Further, the 
eclipse-mapping technique only provides a snapshot of the 
spot-coverage at the instant of the transit. For long-period 
planets, or even stars without transits at all, this technique 
cannot uniquely infer even basic stellar properties, such as 
the rotation period. 

Like so many problems in astrophysics, modelling ro- 
tation modulations due to starspots can be approached 
using either numerical or analytic techniques. Numerical 
techniques are more diverse, allowing one to compute any 
starspot shape, flux profile, limb darkening law, etc one 
wishes. However, they come at the expense of much higher 
computation times. Analytic models are extremely quick to 
execute, often outpacing their numerical counterparts by 
orders-of-magnitude. Such models are challenging to derive 
though and are limited in that they assume fixed proper- 
ties of each spot; for example, their shape is often assumed 
to be circular. However, rotational modulations represent 
a disc-integrated snapshot of a star with unresolved sur- 
face features. For this reason, the size of a spot and the 
flux-contrast are highly correlated since altering either will 
change the amplitude of the resulting rotational modula- 
tions. In this paradigm, modelling elaborate shapes for the 
starspots with dozens of free parameters is unlikely to gen- 
erate a more meaningful model than a simple circular as- 
sumption. Additionally, in the circular model, one should 
interpret the spots as really representing a dark patch or 
cluster of small spots on the surface rather than a perfectly 
circular starspot. 

The diverse range of spots which can be modelled using 
numerical techniques is therefore not a practical advantage 
over analytic models. With this advantage lost, we consider 
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analytic models to be invariably the preferential tool for 
modelling starspots since they will execute with dramati- 
cally quicker computation times. 

1.4 Current Analytic Models 

Analytic models of rotational modulations due to starspots 
have existed for decades. The foundational paper comes from 
lBuddindl| 19771 ) . who describe an analytic model for the rota- 
tional modulation due to multiple non-overlapping circular 
starspots. An alternative derivation of an essentially identi- 
cal model is presented in iDorrenl (1987). 

These models have been successfully applied to numer- 
ous studies of starspots. We highlight the recent work of 
the MOST s pace-telescope in de tecting differe ntial rotation 
on e Eridani (jCroll et al.ll2006bl ) and k 1 Ceti (jWalker et all 
l2007h . In both of these cases, the authors make use of a 
Markov Chain Monte Carlo (MCMC) algorithm to regress 
the data, which yields Bayesian inferences of the parame- 
ter posteriors and correlations. The code, called StarSpotz 
ijCroll et al.ll2006al ). includes parallel tempering to locate the 
global minimum in the inevitably complex parameter land- 
scape. Bayesian inference techniques, such as MCMC, offer 
significant improvements in the statistical interpretation of 
modelling starspots, but come at the cost of being inherently 
computationally expensive. For this reason, analytic models 
are highly prized due to their unmatched computational ef- 
ficiency. 

Despite the successes of the lBuddinj (|l977f ) and lDorrenl 

(1987) model, there are several areas for improvement. 
Firstly and perhaps most critically, the models are limited 
to a linear limb darkening law which is ge nerally a poor 
description stellar specific intensity profiles. IClaretl (|2000h 
remark that the most accurate limb-darkening functions 
are the quadratic and "non-linear" laws, both of which are 
widely used i n the exoplanet comm unity for example. In- 
dee d, recent ly Nutzman et alj l|201ll ) , who also made use of 
the IDorrenl l| 19871 ) model, remarked on how extending the 
model to non-linear laws would be a significant improve- 
ment. Given th e dramat i c incr ease in photometric precision 
since the era of lBuddind (Il977f ) to the Kepler-era,, this point 
is not just pertinent but imperative to address. 

Secondly, there currently exists no partial derivatives 
for the model, which would lead to a further significant im- 
provement in regression analy sis. One of the obst acles in 
achiev ing this goal is that the IBuddind |l977l ) and IDorrenl 
(|l987f ) equation for the model flux is not a single-domain 
function and has multiple cases. This means partial deriva- 
tives would have to be tediously derived for all cases individ- 
ually. Finally, if one sets the goal of deriving partial deriva- 
tives, then it would be advantageous to observers to include 
numerous intrinsic effects in the model, for which their re- 
spective free parameters could also have partial derivatives 
computed. Examples include differential rotation, starspot 
evolution, diluted light and instrumental offsets. 

1.5 Goal of This Work 



For reasons described in £|1.2l fc ill. 31 the principal goal of this 
paper is to provide a revised analytic model for rotational 
modulations due to starspots. This new model has wide ap- 
plications for both stars with and without orbiting planets 



offering nume rous advantages over the IBuddind (|l977T ) and 
IDorrenl (|l987T ) algorithms. Specifically, we highlight the fol- 
lowing key features of the model presented in this paper: 

■ Allows for Ns non-overlapping small starspots, as- 
sumed to be small relative to the stellar radius. 

■ Full non-linear limb darkening of the stellar and spot 
surface is included with vectors c and d respectively. 

■ Differential rotation is included via a latitude- 
dependency including terms in sin 2 $ and sin 4 $. 

■ Starspot evolution permitted using a linear model. 

■ Umbra/penumbra effect may be generated. 

■ M instrumental offsets are allowed for (e.g. quarter-to- 
quarter offsets in Kepler data) 

■ M blended light dilution factors are allowed for (e.g. 
quarter-to-quarter contamination in Kepler data) 

■ Our solution may be expressed as a single-domain ana- 
lytic function. 

■ Consequently, we are able to provide single-domain an- 
alytic expressions for the partial derivatives of the model 
flux with respect to all input parameters, as well as time 
(F'). 

■ We also show how the model may be used to predict 
transit depth variations (TSV) due to non-occulted spots. 

■ We make freely available the new al- 
gorithm in Fortran 90 code, macula (see 
www.cfa.harvard.edu/~dkipping/macula.html). 

In we introduce the model and discuss the various 
definitions and how differential evolution and starspot evo- 
lution is accounted for. In this section, we also provide a 
way to use our model to compute the transit depth varia- 
tions due to unocculted spo ts. In 331 we compare our results 
to that of the lDorrenl (|l987f ) model and show that the small- 
spot approximation used in this work is accurate to within 
-~100 ppm for spots of angular size < 10° . Further, the model 
has a maxi mum error which is a n order-of-mag nitude less 
than that of lBuddind l| 19771) and IDorrenl (|l987V ) for a non- 
linearly limb darkened Sun-like star observed in the Kepler 
bandpass, with spots of angular size below 10°. In i|4] we 
demonstrate an application to a prev iously-studied exam - 
ple, MOST observations of k 1 Ceti bv lWalker et all (|2007f l 
using a multimodal nest ed sampling algorithm, MultiNest 
|Feroz et al.ll2008l . l2009l ). Discussion of the key highlights is 
provided in 



2 THE MODEL 
2.1 Assumptions 

For clarity, we list the assumptions made throughout this 
work below: 

(i) All starspots are circular and lie on the plane of the 
stellar surface 

(ii) Starspots never overlap one another 

(iii) Starspots are small relative to the stellar radius 

(iv) Each starspot is grey and has a uniform temperature 

(v) The star is a sphere with projected circular symmetry 
(i.e. no gravity darkening) 

We do not claim that these are necessarily physically 
true statements. The function of these assumptions is that 
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they are broadly reasonable and allow for a self-consistent 
analytic solution for modelling both in- and out-of-transit 
starspots. 

Out of all our assumptions, the one which is most likely 
to impose practical restrictions on the application of our 
model is the small-spot approximation. One may reason- 
ably question why such an assumption is indeed required. 
By using the small-spot approximation, we may treat the 
surface brightness of the star to be constant under the disk 
of the starspot. This assumption, inspired by the small- 
pl anet approximat i on use d for modelling transit light curves 
in lMandel fc Ago] (|2002f ). leads to dramatically simpler ex- 
pr essions. For example, t he non-linear limb darkening model 
of iMandel fc Agoll (2002) requires hypergeometric functions 
whereas after the authors apply the small-planet approxi- 
mation the most computationally expensive function is arc 
cosine. In the case of a transiting planet, one can see that 
the circular symmetry of the planet is a simpler problem 
than that of the foreshortened starspot, suggesting an ana- 
lytic non-linear limb darkening solution without the small- 
spot approximation would certainly not be computationally 
cheap. 

The small-planet approximation is later shown to be 
dramatically more accurate than the linear limb darkening 
assumption of previous works for almost all feasible spot 
sizes. We estimate a maximum error of 100 ppm for spots of 
angular sizes < 10°, and typically the error is much smaller 
than this. More detailed estimators for the accuracy of our 
small-spot model are provided in Sj3] 



2.2 Definitions 

We define the host star to have Ns starspots on its surface 
labelled by k — 1, 2, ...,Ns — l,Ns- Each spot has a fixed 
angular radius at, which represents the solid-angle of the 
cone swept out from the stellar centre to the stellar surface. 
Further, each starspot has a flux-per-unit-area contrast ra- 
tio, relative to the star, defined by / sp ot,fc. This is essentially 
a proxy for the temperature of the spot and is assumed to be 
uniform within each starspot (but variable between spots). 
Setting /spot.fc > 1 reproduces a bright facula, rather than a 
dark starspot. 

The centre of a starspot has a longitude Afc and latitude 
$fc. These two angles may be combined into the auxiliary 
angle, /3j., defined as 

Pk = cos -1 |^cos /, sin <l?fe + sin /, cos <J>fc cos AfcJ , (1) 

where /* is the inclination of the star. Starspots may 
migrate over time from a reference location due to the stellar 
rotation. We assume that no migration occurs in latitude but 
linear migration is permitted in longitude via 



A fc = A rcfifc + 



$fc = 'I'r 



2ix{t — t re f,fc) 
P*,h ' 



(2) 
(3) 



where P t ,k is the time for the k th spot to undergo a 
change of 2n radians in longitude and t re f,fc is an arbitrary 
reference time when A = Arcf^. We stress here that modify- 
ing our code to include latitude migration is trivial but the 
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Figure 1. An example of our linear starspot evolution model. 
We plot the size of the starspot in units of a max as a function 
of time. The gridlines (left-to-right) mark the end of ingress, the 
instant t max , and the start of egress. 



partial derivatives returned by the algorithm are only valid 
under the above assumption. Due to differential rotation, 
this period varies for each spot and here we assume a simple 
latitude-dependence for the differential rotation of 



P,.fc = 



1 — K2 sin 2 $fc — K4 sin 4 $fc ' 



(4) 



where Peq is the rotation period for the equator of the 
star and K2 and K4 are coefficients of the differential rota- 
tion profile. Additionally, a starspot may evolve via a linear 
growth/decay model of the angular size via 



Qfcfti) 

Q^max.fc 



= Z^ 1 [AtiH(Ati) - At 2 H(At2) 
£: A r 1 [Ai3H(At 3 ) - Ai 4 H(Ai 4 ) 



and using 



. Z/fc 
Aii = u — t max ,fc H — - — hlfc, 



A£2 — ti — tmax.fc + 
At3 — ti imax,fc 



2 ' 
~2~' 



(5) 

(6) 
(7) 
(8) 
(9) 



where a max ,fc is the angular size of the k th spot at a 
reference time t ma x,fc, L k is the "lifetime" of the spot (tech- 
nically the full-width-full-maximum) and X% & £k are the 
ingress & egress durations of the spot's growth profile. H(a;) 
is the Heaviside Theta step-function. An illustrative example 
of our starspot growth/decay model is shown in Figure [T] 

For simplicity, macula defines the reference times t te {^ 
to be equal to i max ,fc, although one may change this def- 
inition without affecting the validity of the returned par- 
tial derivatives. Finally, we employ the f our-coefficient non- 
linear limb darkening law o f lClaretll|2000h . where the specific 
intensity of the star is given by 
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n/2s 



(10) 



where c„ are the limb darkening coefficients, n = 
cos O = Vl — r 2 , ^ r ^ 1 is the normalised radial co- 
ordinate on the disk of the star. We employ the definition 
of a normalised l imb d arkening coefficient, cq, as utilised by 
iMandel fe Agoll (|2002l ) where Co = 1 - c x - c 2 - c 3 - c 4 . 



2.3 Solution 

A detailed derivation of the model presented here is provided 
in Appendices [B] and [C] To summarise, the analytic so- 
lution for the model flux from Ns non-overlapping circular 
starspots may be expressed as 



mod - 2^ Um m I B m F(a = 0, 

m — 1 \ v 



/3) 



, (11) 



where t/ m is the instrumental offset of the m th data set 
(a normalisation factor for each Kepler quarter, for exam- 
ple), B m is a blending factor for each data set (or quarter) 
and n m is a box-car function defined by 



Ilm (J*i , ?start,m , ?end,77i ) = H(t - T st 



H(t — T eIK j,., 



(12) 



In the above, T s tart,m is the start of the m th data set and 
7cnd,m is the end of the m th data set. The F(a, (3) function 
is given by 



and we further define = CiPk — ct k ) and = 
Cifik + afc). We provide some typical examples generated by 
macula using random input parameters for four realisations 
of a 5-spot model in Figure [2] 



2.4 Generating Umbra/Penumbra 

Sunspots can manifest with umbra and penumbra, or sim- 
ply as "naked" spots. Formally, our model assumes naked 
non-overlapping spots. However, our algorithm macula does 
allow one to place spots on top of one another too. Doing 
so allows one to generate umbra/penumbra effects via a su- 
perposition. 

To generate a single starspot with an umbra and penum- 
bra of angular radii a um bra and a pcmlm bra respectively, one 
simply generates two spots of these sizes. If the umbra has 
a flux contrast of / um bra and the penumbra has / pe numbra, 
then the two spots generated will have {a, / ap ot} equal to 

{^p? ./penumbra} and {o?^, /"penumbra /umbra}- 



2.5 Transit Depth Variations from Unocculted 
Spots 

It is well-known that an eclipsing body which occults a 
starspot leaves a s ignificant imprint on the transit profile 
l|Rabus et al.ll2009T ). Non-occulted dark spots also affect the 
transit indirectly via an a mplification of the apparent transit 
depth jCzela et alj |2009); so-called transit depth variations 
(T(5V). This occurs because the planet transits a non-spotty 
region where more flux is concentrated and thus more of the 
total flux is blocked out by the eclipse. The observed transit 
depth is defined as: 



F(a,P) = l-Y i 



n + 4 



4(c n — e2n/ S pot, fc) C_ffc 



n + 4 



C-,fc - C+,fc + ^C+.fc.C- 



(13) 



where S x lV is the Kronecker delta fucntion and c = 
{co, ci, C2, C3, ca} t and d = {do, di, dz, da, di} T describe the 
non-linear limb darkening coefficients of the stellar surface 
and spot surface respectively. The function Ak defines the 
sky-projected area of the k th starspot and is given by: 

A k (ci k ,l3k) = R^cos _1 [cosa fe csc/3fe] 

+ cos/3 fe sina fc H fc - cosa fc sin/3 fc * fc j , (14) 
where we use 



Sfc = sin ctk cos 



cota fc cot/3fc], 



*fc = \J\~- cos 2 a k esc 2 p k . 



(15) 
(16) 



Finally, Equation 1131 includes the ( function, which we 
define as: 



C(as) = cobxH(x)H(^--x) + H(-x), (17) 
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For an unspotted star, this yields: 



lim £ bs = p = 5, 

f.V— >() 



(18) 



(19) 



where p is the ratio of the planet to star radius, Rp/R*. 
To derive the T<5V effect, let us first consider that the spot 
is a bright facula. In this case, the spot actually increases 
the total amount of flux emitted by the star - it provides 
some extra flux Fextra- This extra flux must be given by 



F cxt ra = F(tX, (3) - F(cX = 0,f3). 



(20) 



An observationally equivalent scenario would be to con- 
sider this extra flux as originating from a spatially unre- 
solved background star. This well-known scenario is often 
dubbed a "blend" because the extra flux source is uneclipsed 
and th us the total eclipse depth is diluted due to the blend 
source. Kipping fe Tinettil |2010l ) showed that a blend source 
changes the transit depth via: 



Sobs — — 

B = 



B' 

-P* H~ -Pcxtra 



(21) 
(22) 



The above allows for a simple calculation of the T<5V 
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effect. One additional effect we can include at this point is 
genuine background/foreground blend sources with a blend 
factor Bm- For the bright facula then, the observed transit 
depth becomes 



5ahs 



F(a = 0,/3) 



1 



(23) 



For a bright facula, F(a, (3) > F(a = 0, 0) and thus 
5obs < S (for B = 1) i.e. the transit depth becomes shal- 
lower due to the "blend" source, as expected. This logic is 
easily extended to dark starspots and accurately predicts 
the Tc5V effect, except that dark uncocculted spots cause 
deeper transits. The analogy of a spatially unresolved back- 
ground star becomes unphysical in that the background star 
now emits negative flux, but this is beside the point. One 
key conclusion is that unocculted facula behave as "blends" 
and unocculted spots behave as "anti- blends" . 

Another subtlety is that the above that the effect is 
purely due to blindly normalising the data by the local base- 
line, which is affected by rotational modulation. "Astrophys- 
ical detrending" of a continuous photometric time series, in 
this case using a spot-model to detrend the data, would elim- 
inate any apparent depth variations. However, ground-based 
observers are usually only able to obtain a small amount of 
data either side of a transit event leaving no alternative ex- 
cept to blindly normalise the data. In fact, even space-based 
transit observations are almost always blindly normalised 
using polynomials, running medians or linear trends, macula 
therefore offers the opportunity the astrophysically detrend 
photometry. 

In the typical case of blindly normalised data, Equa- 
tion [221 allows one to fit a set of transit depth variations with 
a spot model using macula. Alternatively, one may wish to 
make causal predictions of the T<5V effect based upon out-of- 
transit rotational modulations. We stress that the equation 
is only valid if the spots are unocculted. These T<5Vs may 
occur in time due to rotational modulation (see examples in 
Fig [5J , or even in wavelength due to the chromatic nature 
of spots. Indeed, correcting for the chromatic T5V effect is 
crucial in ac curate interpretati on of exoplanet transmission 
spectra (e.g. iDesert et al1l2009l ). A more detailed discussion 
of the applications of T<5Vs is presented in N5.3I 

macula directly returns the function (S b s /S) at all 
times, ti, which are inputted. This feature is on/off switch- 
able so that one may either choose to not use the feature or 
perhaps just input one (or a few) value(s) of ti, such as the 
time(s) of transit minimum. Partial derivatives of this func- 
tion are not provided although may be easily computed since 
macula evaluates the partial derivatives of At and F mo d- 



3 COMPARISON TO THE DORREN (1987) 
MODEL 

3.1 Overview 

As discussed earlier in HI. 41 the most wid ely used analytic 
mo dels for starsp ot modelling come from iBuddind (119771 ) 
and |Porrenl(ll987i ). The models are es sentially ident ical and 
so we will refer to comparison to the IDorrenl (|l987l ) model 
only from here-on-in for brevity. 

The main difference between our model and that of 



IDorrenl <|l987h is that our deriv ation as s umes the starspot 
is small i.e. sin a < 0.1, whereas IDorrenl (|l987f ) did not. By 
making this assumption, we have derived a full no n-linear 
limb d arkening treatment for starspots, whereas the IDorrenl 
(1987) model is limited to a simple linear limb darkening 
law only. Further, our model is amenable to the inclusion of 
spot-crossing events due the parametric form of the expres- 
sions describing the arcs along the starspot rim and bulge. 

Since the two models essentially only differ in their 
treatment of limb darkening, one should expect them to be 
exactly equivalent for the case of a uniform brightness star, 
which we easily verified through numerical experiments. 
However, for the case of a limb darkened star, one might 
ask, for what spot size doe s our m odel significantly deviate 
away from that of IDorrenl (|l987l )? We will investigate this 
question in the following subsections. 



3.2 Model Error due to our Small-Spot 
Approximation 

Let us define the model flux, as predicted by this work, as 
F'mod (as used throu ghout). Let us further define the model 
flux as predicted by IDorrenl (119871 ) as F mo d,Dsr- If we as- 
sume a linearly limb d arkened star, t hen the only difference 
between the model of IDorrenl (|l987r ) and that of t his work 
is that we assume a small-spot and IDorrenl (| 19871 ) do not. 
Therefore, direct comparison between the two models for 
linear limb darkened stars yields the error in our model of 
assuming a small-spot. Accordingly, one may write that the 
relative error in our model is: 



AF mod 



I -Fmod,D87 — -Fmod(ci = C3 = C4 = 0; C2 = «L ) I 



Fir 



(24) 



For simplicity, we assume the limb darkening of the spot 
and star are equivalent and set the spot-star contrast, / spo t, 
to be zero (a black spot). Numerically evaluating AF mo d 
over the domain of interest reveals the error is maximised 
when p = a. Therefore, we define AF™ ™ = AF mod (/3 = a). 

The function AF^d grows with both ui and a, tending 
to zero when they both equal zero, as expected. For a Sun- 
like st ar (T eff = 6000 K, log 5 = 4.5 dex, [M/H] = 0). IClaretl 
l|201ll ) estimate that the best fitting linear limb darkening 
coefficient in the Kepler bandpass is ul — 0.5733. One may 
now set AF™od to some desired tolerance level (e.g. the noise 
level of the data) and solve for a i.e. the maximum spot size 
which leads to model errors below the tolerance level. We 
plot this function in Figure [3] (solid-line). 

Since a Mk? p = 12 star has a typical noise of ~ 50ppm 
per long-cadence measurement (note that most Kepler tar- 
gets are fainter than this), we estimate that the maximal 
error of our small-spot approximation is below that of a 
typical Kepler measurement error for starspots of angular 
radius a < 7.6°. This corresponds to a spot coverage of 
< 1.7%. Note that the m odal spot cov erage of stars in the 
Kepler sample is ~ 1% (|Basri et all 1201 ll ). A 1.7% spot 
co verage roug h ly co rresponds to V rns = 0.83 on Figure 4 
of iBasri et all (|201ll ). which can be seen to encompass the 
majority of periodic variables. Given the conservative as- 
sumptions of using a relatively bright 12 th magnitude star 
(most Kepler targets are fainter), the fact that we assume 
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Figure 2. Examples of light curves (solid) and TSVs (dashed) for four randomly generated scenarios using macula. Each simulation 
assumes 5 spots with random properties, including random angular sizes between 0° and 10° and Sun-like non-linear limb darkening. 
Even with 5 spots, the photometric behaviour can be highly complex due to spot evolution, which is also randomly generated in all four 
cases. 



only a single spot is responsible for the entire spot coverage 
(making the spot as large as possible) and the fact that the 
error derived is the maximal error rather than the typical 
error, we conclude that the large majority of spotted stars 
within the Kepler sample are appropriately modelled by the 
expressions in this work. 



3.3 Model Error due to a Linear Limb Darkening 
Law Assumption 

For larger spots, the modelling error becomes larger due to 
our sm all-spot assum ption and thus an observer may opt to 
use the [5 orrc r\ (|l987T l model instead. However, we point out 
that this model assumes a linear limb-darkening law which is 
somewhat unphysical in itself. The question therefore arises, 
at what point is the error in assuming a large spot with linear 
limb darkening better than assuming a small spot with non- 
linear limb darkening? 

We have already calculated the error due to the small- 
spot assumption, assuming the star is perfectly described 
by a linear limb darkening law (A_F mo d). We may similarly 
define an error in assuming a linear limb darkening law when 
the star is really described by a non-linear law: 



A IP mod,D87 mod /otr\ 

^i-Tmod.DSy — — \*0) 

mod 

For a star with the same properties as used in the 
previo us example (T cti = 6000 K, log 3 = 4.5 dex, [M/H] 
= 0), IClaretl (|201lh estimate that the best fitting non- 
linear limb darkening coefficients in the Kepler bandpass 
are {a, c 2 , c 3 , c 4 } = {0.3999,0.4269,-0.0227,-0.0839}. We 
also assume a black spot with the same limb darkening as 
the star, as was done for the previous example. Plotting the 
function A_F mo( j,D87 f° r several realisations of a as a function 
of j9, we find the maximal error occurs at /3 = a/2. Thus we 
define AF™ a £ Dg7 = AF mod , D87 (/3 = a/2). 

In Figure O we plot this function along with AF™^ as 
a function of a for t he Su n-like limb darkening coefficients 
computed bv lClaretl (|2000h . The figure reveals that the error 
in assuming a small-spot is substantially smaller than the 
error in assuming a linear limb darkening law for a ^ 10°, 
as one should expect. As an example, sunspots typically have 
angular sizes < 5° and for a — 5° we find AF™££ — 10 ppm 
whereas AiCod 

D87 = 237 ppm i.e. our model is more than 
an order-of-magnitude more accurate. For spots of size a = 
10° we find AiCS = 162 ppm versus Ai^Sd = 981 ppm. 
We find that the significant error in assuming a linear 
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limb darkening law does not become a better approximation 
than the small-spot model until a > 57.6°, for Sun-like limb 
darkening. After this point, the error in our model rapidly 
tends to infinity and becomes untenable. The exact locations 
of the various turnovers and minima in Figure [3] depend 
upon the limb darkening parameters used. Also, the relia- 
bility of the AF™1 Dg7 function worsens for large a since 
the "truth", assumed to be -F mo d itself starts to become er- 
roneous at high a. 

Neve rtheless , it is c lear that our model is more accurate 
than the iDorrenl (jlQ87h model for even a stars pot equiva- 
lent t o the largest spot ever detected (a ~ 30° ; IStrassmeierl 
Il999l ). Despite this, we would urge observers to use a numer- 
ical approach for such large spots since the model errors are 
significantly greatly than typical measurement errors. Spots 
of size a < 10° should be well-described by the analytic 
model presented in this work. 



4 AN EXAMPLE APPLICATION TO k 1 CETI 

4.1 MOST observations of k 1 Ceti 

k Ceti is a relatively nearby G5 dwarf 9.1 pc from the Solar 
System. The star is notable for having a fairly rapid rota- 
tional period of ~9 days and for being a bright Sun-like star 
at V = 4.84. MOST observations of k 1 Ceti in 2003 revealed 
the presen ce of two starspots w ith rotation periods of 8.9 d 
and 9.3 d (|Rucinski et alj|2004h - However, this single data 
set was insufficient to uniquely determine the latitudes of 
the spots and thus the differential rotation coefficient, K2 
could not be measured (the authors did not consider the 
4 th -order coefficient K4). 

Subsequently, MOST observed k 1 Ceti two more times 
in 2004 and 2005 in order to gather enough data that a 
uniqu e solution c o uld b e inferred. Indeed, this data, reported 
by IWalker et~aH (|2007T ), was argued by the authors to be 
sufficient to locate a single minimum. The authors made use 
of the StarSpotz l|Croll et al.l 12006a ) algorithm to regress 
the data, which in turn employs the Buddin j (| 19771 ) model 
for starspots (note that this is equivalent to the model of 
lDorrenlll987l ). Starspotz locates the global minimum using 
parallel tempering and derives parameter posteriors using 
the Markov Chain Monte Carlo (MCMC) technique. 

The photometry span three data sets, exhibit differen- 
tial rotation and seven spots over three years (ranging from 
a*, = 5.95° to Qfc = 16. 76°) and thus requir ed considerable 
computational effort by I Walker" et all (|2007l ). For these rea- 
sons, the data make for an ideal test of not only our model 
here but for an alternative regressing technique. 



MCMC . A fu ll dis cussion of n e sted sampling is given in 
ISkillind l|2004h and iFeroz et all (120081 ) and for brevity we 
direct those interested to these works. 

Multimodal nested sampling is an implementation of 
the technique to efficiently search parameter space under 
the a ssumption t h at one or m ore modes may exist in the 
data. IFeroz et all l|2008l . l2009h describe multimodal nested 
sampling in detail, in particular in regard to their pub- 
licly available algorithm MultiNest. MultiNest is used 
by the "Hunt for E xomoons with Kepler" (HEK) project 
jKipping et al.ll2012T l to compare the Bayesian evidence of 
a planet versus planet-with-moon regression. We will here 
demonstrate the use of MultiNest with our starspot model 
for the k, 1 Ceti MOST photometry. Currently, MultiNest 
does not make use of the likelihood partial derivatives and 
so the partial derivatives were turned off in our implemen- 
tation of macula. Since we only fit a single model through 
the data, there is no use of the Bayesian evidence here and 
thus we employ the constant efficiency mode of MultiNest 
at a target efficiency of 1% with 4000 live points. 

4.3 Priors 



In ord er to make a fair comparison to the Walker et al.l 
<|2007l ) result, we make the same assumptions as the original 
authors. Accordingly, we assume the same number of spots 
for each data set i.e. 2 spots for 2003, 3 spots for 2004 and 
2 spots for 2005. The spots are assumed to be non-evolving 
over the course of each data set and have a lifetime which en- 
su res they only e x ist wi thin a single data set, as was assumed 
bv lWalker etafl (|2007l ). We also assume f spot ,k = 0.22 for 
all k and that the differential rotation profile is quadratic 
is nature (i.e. we fix = 0). Finally, limb darkening for 
the spot and the star are equivalent and follow a linear law 
governed by ul — 0.6840. Using these assumptions, we have 
the same number o f free parameters (27) as was used by 
IWalker et all |2007l ). 

The 27 parameters are 7 reference longitudes, Ao.fe, 
7 reference latitudes, <I?o,fc, 7 angular radii, ao,k, one 
equatorial rotation period, Peq, one differential rota- 
tion coefficient, k%, one stellar inclination angle, /* and 
three instrumental offset terms, U m . Rather than label 
the offsets by m = 1,2,3, we use m = 2003,2004,2005 
for each year. Since each year has unique starspots, we 
do use the spot labels k — 1,2,3,4,5,6,7 but instead 
use 2003_1, 2003_2, 2004.1, 2004.2, 2004.3, 2005_1&2005_2. 
These labels more clearly identify the spots as sociated with 
each y ear and follow the labelling notation of I Walker et al.l 
(2007). We adopt the uniform priors for all 27 pa rameters 
with the same range as that of I Walker et"ai1 (|2007f l. 



4.2 Multimodal Nested Sampling with MultiNest 

Nested sampling (|Skilling| [2051 ) is a Monte Carlo method 
which puts the calculation of the Bayesian evidence in a 
central role, but also produces posterior inferences as a by- 
product. Nested sampling is generally considerably more ef- 
ficient than MCMC methods for computing the Bayesian 
evidence of a model fit. For exam ple, in cosmological ap- 
plications, Mukhcri ee et al.l (|2006l ) showed that their imple- 
mentation of the method requires a factor of ~ 100 fewer 
posterior evaluations than thermodynamic integration with 



4.4 Results 

The global maximum a-posteriori model fit is presented in 
Figures H \E\ & \E\ for the data sets in 2003, 2004 and 2005 
respectively. Table [1] presents the posteriors of the best fit- 
ting mode a nd compares t hem s ide-by-side with the results 
reported by I Walker et al.1 (|2007l ). 

As revealed in Table [1] the agreement between the de- 
rived system and spot parameters of n 1 Ceti is excellent 
with marginal differences between the estimates. The resid- 
uals of the fits in Figures [4] [5] & [6] closely match those of 
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Figure 3. Left: Comparison of the maximum model error of the [porreri \l98'\) linear limb darkening assumption (dashed) versus the 
small-spot approximation of this work (solid), for a Sun-like star. As expected, over the range of small-spot sizes, the model presented 
in this work is considerably more accurate. Gridlines mark the point at w hich our model is accurare to 50ppm at a = 7.6° . Right: Same 
as left-panel, except we zoom out to a greater x-scale. The \Dorrei\ 1198$) becomes more accurate than the model presented in this work 
for spots larger than 57.6° (marked with gridlines). At this point, the error in both models is 105 mmag and is arguably unusable in 
either case. Note that the location of the minimum in our model error near 50° is sensitive to the limb darkening coefficients used. 



Table 1. Results from fitting the MOST data of k Ceti using 
the macula model presented in this work and the MultiNest al- 
qorithm. Column 2 shows the 68.3% credible range derived by 
Walk er et al\ tSOOj) (taken from column 6 of Table 3 of that 
work). 



Parameter Walker et al. (2007) This Work 



IWalker et all (|2007f ). This therefore shows that macula cou- 
pled with MultiNest is capable of matching the results of 
StarSpotz. 

Some remaining anomalies in the residuals are evident 
and one may be tempted to inp ut more spots to fit these out. 
However, IWalker et alJ (|2007h specifically caution against 
such a process arguing that the anomalies correlate to moon- 
light contamination and other instrumental effects. 

Although it is not the focus of this work to explore the 
inter-parameter correlations and optimal fitting strategies, 
we here briefly comment on this issue. Our fits reveal the 
strongest correlations between a values associated with the 
same data set i.e. the amplitudes of the signal components 
are correlated. We also find that the equatorial period, the 
stellar inclination and the individual latitudes exhibit mu- 
tual correlations, resulting from the uncertainty in the dif- 
ferential rotation determination. 



5 DISCUSSION 
5.1 Performance 

To test the speed of macula, we generated 1000 time stamps 
of a single synthetic input data set for a random choice of 
the star's parameters. In all cases, full non-linear limb dark- 
ening is employed. We generate a single spot with random 
parameters and call macula 100,000 times to evaluate the 
typical execution time. Every call inputted random star and 
spot parameters in order to obtain a reliable average exe- 
cution time. All simulations are run on a single-thread of a 
Intel Core i7 2.9 GHz processor with macula compiled in g95 
using the optimisation flag -03. 

When macula is called, one may instruct the code 
whether to compute the partial derivatives. With derivatives 
turned off, macula requires 0.59 /is per data point. Turning 
derivatives on yields 6.09 fis per data point. Therefore, the 
act of turning on derivatives leads to a slowing down of the 
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Figure 4. Maximum a-posteriori two-spot model fit to the 2003 MOST data of k 1 Ceti using the analytic model presented in this work. 
Regression performed using Multi Nest in conjunction with the 2004 & 2005 data. Residuals to the fit are offset by 0.94- Figure may 
be directly compared to Figure 4 o n Walker et al. I \200j) . where one can see an essentially indistinguishable result. 



code by a factor of ~ 10.3, for a single-spot model. Note 
that these times include the small overhead of generating 
random system parameters too. 

Increasing the spot-number, we find that the no- 
derivatives call scales linearly with Ns- However, the deriva- 
tives call exhibits super-linear, yet sub-quadratic, scaling of 
Ns' 74 , or roughly N 7 J A . For this scaling, doubling the num- 
ber of spots increases the CPU time by a factor of 3.34. 

Due to its analytic nature, macula performs significantly 
faster than nu merical codes made available in the literature. 
For example, iBoisse et al] (|2012h presented their numeri- 
cal algorithm SOAP and report that generating 10,000 time 
stamps of a single synthetic starspot requires less than 40 s 
(but presumably close to this value). This indicates SOAP 
requires ~4 ms per data point, compared to macula which re- 
quires 0.6 ^s per data point i.e. macula is around 6800 times 
faster than SOAP. There are several points which make a di- 
rect comparison somewhat unfair though, macula does not 
compute radial velocity variations, whereas SOAP does (al- 
though £15.31 shows how radial velocity variations are easily 
generated from the outputs of macula) . Further, the authors 
used a slower 2.33 GHz Intel Core Duo processor for their 
benchmark tests. Nevertheless, it is clear that the difference 



in computation speeds is three orders-of-magnitude, making 
macula a powerful tool in inverse-problems. 

5.2 Benefits 

The analytic model for starspots presented here has several 
advantages which we list here: 

■ An analytic algorithm for modelling photometric rota- 
tional modulation due to multiple circular, grey starspots, 
performing three orders-of-magnitude faster than compara- 
ble numerical codes. 

■ Reproduces light curves with a maximum model er- 
ror an o r der-o f-ma gnitude l ess th an that of the previous 
iBuddind <| 19771 ) and lDorrenl (|l987f ) for a Sun-like non-linear 
limb darkened star observed in the Kepler bandpass, for 
spots of angular size < 10°. 

■ Model accounts for spot contrast, non-linear limb dark- 
ening, differential rotation and starspot evolution. 

■ Includes baseline normalisation parameters for M data 
sets, as well as M blended light dilution factors to aid in 
Kepler analysis. 

■ Computes transit depth variations (T<5Vs) due to un- 
occulted spots. 
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Figure 5. Maximum a-posteriori two-spot model fit to the 2004 MOST data of k 1 Ceti using the analytic model presented in this work. 
Regression performed using Multi Nest in conjunction with the 2003 & 2005 data. Residuals to the fit are offset by 0.97. Figure may 
be directly compared to Figure 5 o n Walker et al. I \200j) . where one can see an essentially indistinguishable result. 



m Partial derivatives of the model flux is provided with 
respect to all model parameters and time, and may be turned 
on/off as desired (see Appendices IDBJEI for derivations). 

■ Code is freely available as a Fortran routine, macula, 
located at www.cfa.harvard.edu/~dkipping/macula.html. 



5.3 Potential Applications 

5.3.1 Rotational Modulation Measurements 

We foresee several possible applications of macula. Firstly, 
measuring the rotational modulation of variable stars may 
be used to determine the rotation period, which may in 
turn constrain the ages of stars with gyrochronology (|Barnesl 
2009). Stars monitored with high signal-to-noise continu- 
ous photometry, such as that from Kepler, may also re- 
veal differential rotation and the stellar inclination angle. 
An example of thi s type of regression is the analysis of 
Walker et al. I i|2007l ) for k 1 Ceti, where the inclination angle 
derived from rotational modulation alone and an analytic 
model for starspots yields a result fully consistent with the 
spectroscopic Vsin/* value. Note that we also reproduce 
this result using macula in this work. In addition, it may 



be possible to measure starspot evolution using the linear 
model employed by macula. 

Rotational modulation analyses using macula are not 
limited to cool stars, which have been most commonly ob- 
served to exhibit such behaviour. There also exists evidence 
for spots on hot stars, such as th e rapidly rotating B star 
HD 174648 l|Degroote et al.ll201lh . Indeed, macula will also 
be applicable for bright spots on hot mass i ve sta rs, such as 
those proposed bv lCantiello fc Braithwaitel (|201ll ) to explain 
observations of late O-type and early B-type stars made by 
CoRoT. 

macula also produces predictions for the transit depth 
variations (T<5V) at any time stamp inputted. This may per- 
mit for the determination of rotational periods from T<5Vs 
alone; highly useful for ground-based observations lacking 
the continuous photometry of space-based observatories. It 
may also be useful in testing whether observed T5Vs are 
consistent with starspots versus som e other hypothesis e.g . 
planetary oblateness with precession l|Carter fc Winnll2010l ) . 

5.3.2 Astrophysical Detrending 

Due to the analytic nature of macula, the code is quick to ex- 
ecute and therefore many find uses in astrophysical detrend- 
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Figure 6. Maximum a-posteriori two-spot model fit to the 2005 MOST data of k 1 Ceti using the analytic model presented in this work. 
Regression performed using Multi Nest in conjunction with the 2003 & 2004 data. Residuals to the fit are offset by 0.97. Figure may 
be directly compared to Figure 6 o n Walker et al. I \200j) . where one can see an essentially indistinguishable result. 



ing of photometry. For example, the PDC-MAP algorithm 
of Kepler is designed to remove instrumental trends but pre- 
serve the astrophysical signal, such as rotational modulation 
due to starspots. In performing a search for transits, or a 
detailed modelling of a known transit, detrending this rota- 
tional modulation is required. Whilst polynomial models or 
harmonic filtering may be used, an astrophysically-grounded 
model, such as macula, offers a viable alternative due to its 
fast execution time. 



5.3.3 Radial Velocity Variations due to Starspots 

We briefly remark that macula may be used to predict radial 
velocity var iations due to stars pots via the F F' method 
described in lAigrain et all (|2012t) . Here, the authors propose 
that radial velocity variations can be reliably predicted from 
flux variations (F) alone. Specifically, the authors argue that 
the flux multiplied by its derivative in time reveals the radial 
velocity variations, macula returns both -Fmod and <9F mod /<9t 
for all input times (derivation provided in Appendix [FJ. 



5.3.4 Correcting Transmission Spectra 

If a planet transits across a star, the atmosphere of the 
planet can reveal chromatic variations in the transit depth 
due to molecular absorption. In this way, transit measure- 
ments can reveal the "transmission spectrum" of an exo- 
planet. If the host star has unocculted starspots, one would 
expect to see chromatic variations in the depths purely from 
spots too, introducing a confounding signal. Further, transit 
depth measurements are often scattered both in wavelength 
and in time meaning that rotational modulation can also in- 
troduce spurious depth variations. Correcting for starspots 
is therefore a major challenge in studying the atmos pheres of 
extrasolar planets. For example. iDesert et alj (|2009T > found it 
necessary to use rotational modulation data of HD 189733 in 
order to correct Spitzer measurements of the planet's transit 
depth. 

macula offers a self-consistent way of modelling such 
corrections. Observations of rotational modulation may be 
used to directly infer T5Vs, provided the data are in the 
same bandpass as that used for the transit measurements. If 
the bandpasses differ (which is practically speaking likely), 
then one can estimate the depth variations by assuming a 
model for the spectral radiance of the starspots (e.g. black- 
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body). Nevertheless, we point out that obtaining several 
spectra within a few rotation periods of the intended transit 
measurement would be the most ideal way to correct for such 
activity, obviating the need for spectral radiance modelling. 

5.3.5 Exomoon False- Positive Vetting 

Whilst we leave the issue of modelling planet-spot cross- 
ings to future work, it is worth noting that such events may 
resemble exomoon mutual transits and are anticipated to 
be a source of false-positives in the "Hunt for Ex omoon 
with Kepler" (HEK) project (|Kipping et al.1 \201$ ). Even 
without detailed spot-crossing models, macula offers some 
simple tests to compare these two competing hypotheses. 
Firstly, the starspot coverage can be estimated from the 
out-of-transit variability, allowing one to gauge the feasi- 
bility of an observed anomaly being a spot-crossing event. 
Secondly, the derived rotation period from the rotational 
modulations can be used to check whether the light curve 
anomali es are consistent or i n consis tent with such a period. 
Finally, ISanchis-Oieda et al.1 (|2012l ) (see Fi gure 1) have re- 
cently shown that the phase of a transit mid-time with re- 
spect to the rotational modulations ($tra) is related to the 
phase of a starspot with respect to the transit mid-time 
(<l?spot). Any variations not matching this relationship would 
be difficult to explain as being due to a starspot. 
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APPENDIX A: GEOMETRY OF THE SPOT 



Al Position of the Spot 

The basic unit of our model is a circle, which represents a 
starspot, on the canvas of a star's surface. We begin by con- 
sidering a single starspot and show later how the result is 
generalised to multiple spots. We will assume that the star 
is perfectly spherical in what follows. The geometry of the 
spot is characterised by a position and a size. The position, 
which we define as the position of the starspot's centre rel- 
ative to the centre of the star, must be a two-dimensional 
vector given that a surface has a two-dimensional topology. 
An appropriate positional vector would be longitude (A) and 
latitude ($). We define these terms to exist in the range 
— 7r < A < 7r and — 7r/2 < $ < tt/2. 

We initially consider the centre of the spot to be located 
in a Cartesian frame at a location given by the unit vector 
k = {0, 0, 1} T (where we adopt units of the stellar radius). 
In all frames of reference, we consider the observer to be 
located along the z-axis at z — +00. 

The centre of the spot can be described at any longi- 
tude and/or latitude by multiplying the unit vector k by two 
rotation matrices, accounting for longitude and latitude. At 
this stage, we denote the longitude and latitude using the 
notation A and $ respectively, which we dub "apparent lon- 
gitude" and "apparent latitude". This is done in order to 
reserve the symbols A and $ (the true longitude and lati- 
tude) for later when we will account for stellar inclination as 
well. The rotation matrices for apparent longitude and lati- 
tude are defined by the notation M A then M§ , respectively. 

Consider that the action of these two rotation matrices 
leads to a position for the centre of the spot defined by the 
vector Rcentre- Due to the non-commutative nature of linear 
algebra, the order in which one chooses to perform these 
rotations will affect the results. Here we follow historical 
precedent and define: 



Rcc 



M A M ti ,k, 



(Al) 
(A2) 



^centre 




sin A cos <E> 


^centre 




sin $ 


^centre 




cos A cos <!> 



we have 


\x'- ' 

'"rim 




sin a cos ip 


r>' 

■"rim 


/ 

2/rim 




sin a sin ip 




/ 

_ rim_ 




cos a 



Mi 



M = 



cos A 
1 
— sin A 



1 

cos $ 
— sin <l 



sin A 


cos A 





sin $ 
cos $ 



(A3) 



(A4) 



Rcentre = J/centre = sin $ . (A6) 



A2 Size of the Spot 

We wish to define the radius of the spot in terms of a solid 
angle swept out from the centre of the star. Let us define this 
solid angle to be given by a. A spot of solid angle radius tt/2 
radians would reach from pole-to-pole and thus we define 
< a < tt/2. Later when we account for limb darkening 
effects ( i]C2|) . we show that it is necessary to assume < 
q < 7r/4 and this should be interpreted as the hard- limit of 
our model, macula. 

For a starspot with a position vector described by 
Rcentre = k, it is trivial to show that the apparent radius of 
the spot would be sin a. In this frame, the spot appears as 
a perfect circle on the X-Y plane. 



A3 Rim of the Spot 

We define the rim of the spot to be those points which 
lie along the two-dimensional projected perimeter of the 
starspot, when viewed along a vector normal to the stel- 
lar surface and passing through centre of the starspot (i.e. 
when Rcentre = k) . After applying the rotation matrices to 
account a starspot's apparent longitude and/or latitude, the 
position vectors describing the loci of points along the rim 
are transformed too. 

Let us define the position vector of the loci of points 
along the rim, when viewed in the frame such that Rcentre = 
k, by the vector R rim . After accounting for the spot's ap- 
parent longitude and latitude, we use the vector Rrtm. 

For R r i m , the loci of points may be described using para- 
metric equations, taking advantage of the fact the projection 
of the spot is a perfect circle (as described in the previous 
subsection). 



(A7) 



where < p < 2tt traces the loci of all points along the 
starspot rim. We may now apply the rotation matrix M A 5 
to find the parametric expressions describing the rim for any 
apparent longitude or latitude, thereby accounting for the 
fore-shortening effect. 



One may combine the two matrices into a general trans- 
formation matrix, M A ^, given by 



Rrim — -M-A |>^rin 

which may be shown to yield 



(A8) 



cos A — sin A sin $ sin A cos $ 

cos i> sin l> 

— sin A — cos A sin <E> cos A cos <E> 



We use this matrix to determine 



(A5) 



Xii m — sin a cos A cos ip + sin A(cos a cos <3? — sin a sin <E> sin ip) , 

(A9) 

2/rim = sin a cos $ sin ip + cos a sin <E>, (A10) 

Zrim = cos a cos A cos $ — sin tv(sin A cos ip + cos A sin $ sin tp) . 

(All) 
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A4 Bulge of the Spot 

Consider again the frame in which one views the spot down 
the vector normal to the stellar surface and passing through 
the spot's centre. In the model described in this work, the 
spot lives in three-dimensions in a Cartesian framework. No- 
tably, the spot exhibits a bulge due to the curvature of the 
stellar surface. From the perspective of the star's centre, the 
loci of the points on this bulge can be described by two an- 
gles; a radial angle, 8, and an azimuthal angle, v. We may 
define these loci by again starting from the frame in which 
Rcentro = k, and applying rotation matrices appropriately. 
In this simple frame, we define the position vector for the 
loci of the points existing on the bulge as 



Rv 



bulge 



= M„M fl k, 



where we have 



M„ = 



1 
cos 8 sin 8 
— sin 8 cos 8 

cos v sin v 
— sin v cos v 
1 



(A12) 



(A13) 



(A14) 



Here the radial angle, 8, is bound to be —a < 8 ^ a 
i.e. it cannot subtend an angle greater than the solid angle 
radius of the spot. The azimuthal angle has the freedom to 
be —tv < v < n. We use these matrices to determine: 



R-hi 



IV 

■^bulgc 




sin sin v 


/ 

?/bulgc 




sin 6 cos v 


/ 

_^bulgc_ 




cos 6 



(A15) 



Note the dash, which (as before) is used to denote that 
this is derived in the frame not accounting for a spot's appar- 
ent longitude and/or latitude. As was done earlier, we may 
now apply the longitude-latitude rotation matrix, j=, to 
account for any orientation desired: 



R 



bulge 



(A16) 



which gives 



a^buigc = sin A cos $ cos 8 + sin #(cos A sin v — sin A sin $ cos v) , 

(A17) 

2/buigc = cos <E> sin 8 cos v + sin $ cos 8, (A18) 

Zbuigc = cos A cos $ cos 8 — sin 8 (sin A sin v + cos A sin $ cos v) . 

(A19) 

A5 A Useful Simplification 

Due to the circular symmetry of the problem, it is actually 
degenerate to use two angles to describe the position of the 
spot. All that matters is how close to the spot is to the limb, 
regardless as to the combination of longitude and latitude 
responsible. For this reason, we may define any combination 
of these terms using a single " auxiliary ang le" we dub /3 
(in- keeping with the notation of iDorrerJI 19871 ) . 



The angle of interest is the angle subtended between 
the vector k (pointing towards the observer) and the vector 
describing the position of the spot's centre relative to the 
centre of the star R ccntro . Let us define this as the auxiliary 
angle /3. f3 can be found by using the dot-product rule of 
these two relevant vectors: 



Rc 



k= IRc 



|k| cos/3 



(A20) 



Since k is a unit-vector in the Z-direction, then this dot- 
product simply extracts the Z-component of R ccntro . There- 
fore we have: 



/3 = COS 1 [Zcentre] , 

which may be evaluated here to be 
/3 = cos - 1 [cos A cos <I>] . 



(A21) 



(A22) 



/3 may also be thought of as being like a net longitude 
shift at zero latitude i.e. <3? — > and A — > f3. 

Due to the mirror symmetry of the problem, we only 
need consider < /? < tt to derive all possible scenarios. 
The vectors of interest now become, without any loss of 
generality, 



Rc 



■Rrin 



^centre 




sin ft 


D centre 







•^centre 




cos ft 



(A23) 







sin a cos ft cos ip + cos a sin (3 




2/rim 




sin a sin if 


, (A24) 


•2rim 




cos a cos j3 — sin a sin f3 cos ip 







3^bulgc 




I^bulgc 


2/bulgc 






■^bulgc 





sin P cos 8 + cos /3 sin 8 sin v 

cos v sin 8 
cos 8 cos f3 — sin 8 sin /3 sin v 



(A25) 



APPENDIX B: FOUR CASES 
Bl Overview 

Bl.l Case I 

In order to compute the flux from a starspot, we need to 
compute the projected area in the X-Y plane. It can be 
easily seen that four distinct cases exist for the geometry of 
the rim and bulge. The most obvious case is the dominant 
source of flux variations since the spot is nearly face-on. 
For f3 between and some angle close to the limb of the 
star, the loci of points on the bulge lie fully inside the rim 
of the starspot, as seen in the projected X-Y plane. This 
case is trivial to model and the rim expressions may be used 
alone to compute the area of the starspot. Case I is valid for 
< f3 < /Scrit where we are yet to define /? C rit but it can be 
understood to be angle close to the limb of the star. 

B1.2 Case II 

Case II occurs as /3 approaches n/2 from 0. It is defined as 
when the loci of points on the bulge are no longer contained 
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within the projected rim of the starspot. Since the bulge 
has a Z-component, as we rotate round in longitude, this 
Z-component will be transferred into an ever-increasing X- 
component. Eventually, this X-component exceeds the rim's 
maximal X-value at which point "the bulge pokes out of the 
rim". Case II is valid for B CI it < ft < 7r/2. 



B1.3 Case III 

Case III occurs as (3 increases beyond tt/2 i.e. the centre 
of the spot is behind the star. However, a portion of the 
starspot is still in view and causes a flux decrement. It can 
be easily understood that once B exceeds tt/2 + a then the 
spot has fully disappeared behind the back of the star. Thus, 
case III is valid for tt/2 < B < tt/2 + a. 



B1.4 Case IV 

Case IV is simply the case of the spot fully behind the star 
and thus there is no contribution to the model flux. This is 
valid for tt/2 + a < B < tv (recalling that B is defined only 
within the range < 8 < n due to the mirror symmetry of 
the problem). 



B2 Case II 

B2.1 Optimal Bulge Curve 

We here devote a section to case II alone, due to the non- 
trivial nature of solving for its parametric equations. We first 
start by defining the "optimal bulge curve" . For any Y co- 
ordinate of a point lying within the bulge, the optimal bulge 
curve is the corresponding X co-ordinate which maximises 
X. It is the curve which seems to extend furthest to the limb 
of the star, as seen in the transformed frame. Since all loci 
on the 2D surface of the bulge are defined by two paramet- 
ric terms (9 and u), it should be clear that the parametric 
equation of the optimal bulge curve will require only one 
term; either 9 or v, but not both. We arbitrarily choose here 
to define our optimal bulge curve purely in terms of v. 

The optimal bulge curve also exhibits the greatest sep- 
aration from {X, Y} = {0,0}, relative to all other loci on 
the bulge. Thus we expect that Xbuige + 2/buige is maximised 
and so: 



cK^bulgc + J/bulgc) 



0. 



(Bl) 



Solving the above for cos 9 yields two solutions, only 
one of which is the maximum: 



COS ^optimal — 



2 sin 8 sin v 



x/3 + cos 2/3-2 cos 2v sin 2 



(B2) 



where we only consider the range < v < tt/2 and 
< B < tt/2 here (the latter due to the case II conditions 
and the former due to symmetry about the X-axis). This 
yields the following parametric expression in the X-Y plane: 



3-optir: 



2/optir 



2 sin V 



^3 + cos 2/3-2 cos 2v sin 2 B ' 

2 cos v cos B 
y/3 + cos 2/3-2 cos 2v sin 2 8 ' 



(B3) 
(B4) 



Evaluating the equation for x op timum at v = reveals 
^optimum = 0. Thus, when v = the optimal bulge curve 
intersects the a;-axis, although we note that at this point 
the corresponding ^optimum point may be exceed a and thus 
may not truly exist on the bulge. However, it reveals that x 
increases as v increases form to 7r/2. 



B2.2 Intersection of Optimal Bulge Curve and the Rim 

For case II, where B cr it < /3 < tt/2, there exists a certain 
point where the optimal bulge curve intersects the starspot 
rim. We denote this location as {^intersection, ^/intersection}- 
The location corresponds to a unique parametric loca- 
tion along the rim, (^intersection- Similarly, there exists a 
unique parametric location along the optimal bulge curve, 

^intersection • 

Let us deal with (^intersection first. Since the optimal 
bulge curve extends outside the rim, this location can be 
shown to occur when the rim's X-Y distance from the ori- 
gin is maximised i.e. when a; 2 m + i/ 2 im is maximised. We 
therefore must solve the following expression for v. 



^(-^rim ~j~ 2/rim) 

dip 

which may be shown to yield: 



COS (^intersection = COt O COt /3 



(B5) 



(B6) 



The intersection point along the optimal bulge curve 
can be found by minimising the distance on the X-Y plane 
between the optimal bulge curve and the rim. Therefore, we 
must solve the following expression for v. 

~q ^[^optimal Xrirni^P = (^intersection)] 

+ [^optimal — 2/rim(</3 = (^intersection)] 2 ^ = 0, (B7) 

which yields the following solution: 



COS 2 ^intersection = 1 — COt 2 Ct COt 2 ft. 



(B8) 



Feeding this back into the expressions for the optimal 
bulge curve, we locate the Cartesian co-ordinates of the in- 
tersection point: 

^intersection = COS Q CSC B, (B9) 
^intersection = sin Ci\J 1 - COt 2 Q COt 2 8, (B10) 

where it is again understood this is for the range < 
v < tv/2 only. 
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B2.3 The Critical Angle, /3 CTit 

As discussed earlier, ^optimal = for v = and increases 
up to a maximum at v = tt/2. Similarly, by definition the 
parametric expression for :r r im is maximised for <p — 0. Case 
II is only valid for a bulge which pokes out of the rim and 
its boundary will occur for a; op timai(^ = ir/2) = x rirn (<p — 0). 
Solving for /3, we find: 

cos^ cr it = sin a, 

/3 C rit = tt/2 - a. (Bll) 

This therefore proves an intuitive point. The optimal 
bulge pokes out of the rim when the rim hits the edge of 
the star. At this instant, the starspot rim makes contact 
with the projected rim of the star and the optimal bulge is 
just a single point at {x,y} = {1,0}. As /3 becomes larger, 
the starspot rim gradually disappears behind the back of 
the star and the optimal bulge curve spreads out along the 
projected rim of the star. 

B2.4 Projected Area of the Starspot: I. The Rim 

The projected area of the starspot, for case II, can be 
thought of as the sum of the projected area bound by the 
rim and that of the bulge poking out, with the transition 
occurring at the intersection points derived. The rim there- 
fore bounds an area between (/^intersection < f> < (27r — 
^intersection ). We consider here the area above the a;-axis 
only, which can later be simply doubled due to symmetry 
about the x-axis. 

We start by re- writing the expression for x T im(<p) to 
make <p the subject: 

<fi(Xrim) = COS _1 ^ CSC Of SCC /3(:Trim — COS Ct SU1 /3)^ (B12) 

We may now replace the ip in j/rim(v) to obtain 

2/rim (^rim) : 



2/rim (-^rim) 

= sina^/l - ( esc a sec /3 — cot a tan /3) 2 

(B13) 

The area bounded by the rim is therefore given by: 



/"^intersection 
^rim — 2 / ^/rim (^rim) d^rim , 

J Zrim (•/>=?■") 

1 . 2 

Arim = —J= COSQ\/- COS 2a — COS 2/3 cot j3 



(B14) 



+ 



cos 2 /3 sin 2 a + sin 1 |^ cot a cot /jj cos /3 sin 2 a. 

(B15) 



B2. 5 Projected Area of the Starspot: II. The Bulge 

We now need to repeat this process for the optimal bulge 
curve. One may re- write the expression for a; op ti m ai(i / ) mak- 
ing v the subject: 



^(Xoptimal) = COS 



'optimal 



2 ^'optimal ~h ^optimal COS 2/3 



Feeding this into the expression for j/ ptimai (^) m order 
to obtain y op timai (^optimal) we obtain the simple solution: 



^/optimal 



1 — r 2 

A ^optimal" 



(BIT) 



Once again, this result proves the same inuitive result 
we saw earlier. Specifically, the optimal bulge curve lies along 
the projected rim of the star itself. The bounded area is given 
by: 



-<4onti 



optimal 



/ ^optimal (^V 2 ) 
^optimal (^optimal) d^optimalj 
- intersection 

(B18) 



^optimal = COS [cOS « CSC 

cos a 



\/2cot 2 /3^-cos 2a -cos 2(3 



+ 2 sin a tan (3\/ — cot 2 a cot 2 (3 + cos 2 (3 esc 2 a 



(B19) 



B2.6 Projected Area of the Starspot: Total 



Combining these two results together, we obtain the area of 
a circular starspot under case II conditions: 

An(a,P) = cos _1 [cos a esc (3] 



+ sin a 



cos/3sina(7r — cos 1 [cot a cot /?] ) (B20) 



(B16) 



— cos a tan {3\[— cot 2 a cot 2 (3 + cos 2 (3 esc 2 a 



(B21) 

B3 Case III 

B3. 1 Edge Bulge Curve 

For case III, we have 7r/2 < (3 < 7r/2 + a. Here, the centre 
of the spot is out-of-view, hidden behind the star. Despite 
this, a portion of the spot's surface remains at z > and 
thus is still visible. For case II, we defined an optimal bulge 
curve which tracked the curve of interest. Similarly, here we 
define the "edge bulge curve" to the perimeter of the bulge 
still in view when case III conditions remain in effect. 

The edge bulge curve is much easier to define that the 
optimal bulge curve. For the range tt/2 < (3 < n/2 + a, 
it is simply given by maximizing 9. Since 9 is bound to be 
< 9 < a, then # e d g c = a. Thus, the parametric equations 
describing the edge bulge curve are: 

■£edge ■Ebulge(^ ^0) 

= cos a sin (3 + cos (3 sin a sin v, (B22) 

2/cdgc = J/bulge(6> = a), 

= cos v sin a, (B23) 
= Zbuigc(# = a), 

— cos a cos (3 — sin a sin (3 sin v. (B24) 
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B3. 2 Boundary of the Edge Bulge Curve 

The edge bulge curve intersects the stellar rim when the 
quadrature sum of the X and Y components equals unity. 
Therefore, we may find the v value of this location, which 
we dub boundary, by solving the following expression for v: 



2 4- 2 — 1 

■^edge "T" Vcdgc 



which yields: 



2 cos 2a + cos 2/3 

COS ^boundary 



(B25) 



(B26) 



2 sin 2 a sin 2 f3 

Plugging the above into our expressions for R c d gc yields 

R b 

oundary ■ 



^boundary = COS OL CSC /3, 

CSC P I 

J/boundary = l=r V ~ COS 2a - COS 2fj, 

v2 

^boundary — 0. 



(B27) 
(B28) 
(B29) 



The right-most X-point occurs when we cross the X- 
axis i.e. when y c d g c = 0. It is trivial to show this occurs for 
v — tt/2 and correspondingly a; c dgc(^ = 7r/2) = sin(/3 + a). 

B3.3 Area Bounded by the Edge Bulge Curve 

Taking the expression for s c dgc(j / ), we may re- write this to 
make v the subject via: 

//(xcdgc) = sin -1 |^csc q sec /3(a; c dgc — cos a tan ,9)J . (B30) 
We may feed this into y c dge(i y ) to obtain y e dge(£edge): 



J/edget^edgo) = sin q \J 1 — (x csc ct sec P — cotatan/3) 2 . 

(B31) 

Due to the concave nature of the edge bulge curve, the 
area of the loci of points on the bulge only is defined by: 



. 1 / , - | / sj\~x> dx 

'^boundary 
f ZcdgcO^/ 2 ) \ 

2 1 /' ycd gc (a;cdgc)da; cdgc . (B32) 

^boundary / 

Finally, one may express this purely as a function of a 
and /3: 



Alll(a, yS) = SB) C"SI! CM' 



1 |^ cos a csc /?j 
+ cos P sin 2 q cos" 1 |^ — cot a cot /jj 
— cos a sin /3 \J\ — cos 2 a csc 2 ft (B33) 

B4 Cases I & IV 

B4.1 Case I 

For completion, we here briefly derive the expressions for 
cases I and IV. Case I has the spot fully in view at some angle 



/3 where < (3 < /3 C rit- The relevant parametric equations 
are the rim expressions derived earlier. The area may be 
found to be: 



Ai(a, P) — ty sin acos/3 



.2 Case IV 



im \ j 
J/rim— 7: I <i(f 



dtp 



dtp 



(B34) 



Case IV is for (-k/2) + a < P < -k and corresponds to the 
spot fully out-of-view behind the star. The case trivially has 



A IV {a,P) = (B35) 



APPENDIX C: MODELLING THE LIGHT 
CURVE 

CI For a Uniform Brightness Star 

Cl.l Generalising to a Single Domain Function 

It can be easily shown that these ex pressions prod uce the 
same light curve profile predicting by iDorrenl f|l987| s | in the 
absence of limb darkening and a black spot. 

The expressions for An and Am possess some simi- 
larities in form and are of course continuous at the point 
P — n/2. This led us to investigate if the two equations are 
equivalent to some simplified form. We found the following 
expression describes both An and Am: 

A(a, P) — cos -1 |^cos«csc/3j 

+ cos P sin aH — cos a sin P^, (CI) 

where 



3 = sin a cos 1 [— cot a cot P] , $ = yl - cos 2 a csc 2 p. 

(C2) 

Encouraged by this, we tried plotting the function in 
the range < P < /3 cr it- However, A becomes complex in 
this range. We therefore only considered the real part. It 
is easy to see by example that the real part of A perfectly 
maps the Ai function. 

A final success of A comes from considering the case IV 
range i.e. (it/2) + a < P < ir. Here the real part of A goes 
to zero but the imaginary component gradually increases. 
Thus, by plotting the real part of A only, we can reproduce 
all four cases with a single function across the full domain 
of < P < n. Thus, we have: 



A(a,P) =R[A]. 



(C3) 



The advantage of using this single-domain function is 
that we can define an analytic Jacobian and Hessian matri- 
ces, which are useful in expediting regression of photometric 
data, macula will provide the Jacobian, but not the Hessian 
to save computation time (although it may be extended to 
perform this function too due its analytic form). 
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CI. 2 Model Flux 



For a single-rotating spot on a uniform brightness star, the 
flux from a star (F) can be computed using: 



F{a, /?) = (tt - A)T, + AT spot , 



(C4) 



where J- denotes flux-per-unit-area and A is the area of 
the spot. Photometric observations are usually normalised 
to some arbitrary value. A suitable choice here is the flux 
from the star in the absence of a starspot i.e. F(a = 0, /3). 



F, 



gKg) 

F(a = 0,/3) : 



l--(l-/spot) 
7T 



(C5) 



where / spot = Fspot/F* and is the flux-contrast of the 
spot relative to the star. It may also be thought of as a proxy 
for the temperature of the spot. For Ns non-overlapping 
starspots labelled k = 1,2,..., Ns — l,Ns, this can be ex- 
tended to: 



-Fmod(a,/3) = 



F(a = 0,(3)' 



(C6) 



which yields 



F mod (a,f3) = l--V4(l-/ spotii ). (C7) 

7T z ^ 



41 

Ispot(r) = /spot f 1 - d n (l 



/' 



re/2- 



(C9) 



Note how we assume / spo t is not a function of position 
on the star's surface or equivalently the spot's radial angle, 
9. In other words each spot has a uniform temperature, al- 
though the temperature may vary between spots. Also, £12.41 
describes how umbra/penumbra may be generated using our 
model allowing for a more complex profile. 

For a single rotating spot, we will denote the model flux 
as being composed by the following components: 



F(a,/3)=Ftot^-F o] 



+ F S 



pot ■ 



(CIO) 



The obscured and total flux components are computed 
as one would do so for a transiting planet model. In this sce- 
nario, there are four principal cases, as shown in Table EU 
The table requires we define an angle at which point the spot 
no longer covers the centre of the star (M3-M9 boundary). 
When does this occur? 

In order to simplify the problem, let us assume it is not 
possible for a spot to both cover the stellar centre and to 
exceed the angle /3 C rit- This is equivalent to assuming < 
a < 7r/4. For the loci of points along the rim of the starspot, 
the locus which is closest to the sky-projected stellar centre 
has a position {x,y} — {sin(/3 — a),0}. Therefore, when 
sin(/3 — a) > 0, the spot no longer covers the stellar centre. 
The spot therefore no longer covers the stellar centre once 
P > a and this is the critical angle of interest required for 
deriving our limb darkening model. 



C2 For a Limb Darkened Star 

C2.1 The Mandel-Agol Gases 

In this work, we will assume that the size of the spot is small 
relative to the size of the star. This approximation allows us 
to easily write down analytic f unctions for the l ight c urve 
and follows on from the work o f lMandel fc Agoll <|2002h and 
iKippin el (|201ll ). Spe cifically, we will use t he small-planet 
approximation from iMandel fc Agoll {2002) and utilise 4- 
coefficient non-linear limb darkening. The flux from a star 
in the absence of starspots is therefore modelled via Equa- 
tion [TD] provided earlier: 

4 

J,(r) = l-y>„(l-M" /2 ), (C8) 

n=l 

where c„ are the limb darkening coefficients, /i = 
cos G = yl — r' 2 , ^ r ^ 1 is the normalised radial co- 
ordinate on the disk of the star and I*(r) is the specific 
intensity as a function of r, with 1,(0) = 1. 

Limb darkening is present over the entire viewable sur- 
face of the star, including those portions which are covered 
in starspots. However, due to the different temperature and 
opacity of this surface, the limb darkening coefficients can- 
not be assumed to be necessarily the same as that for the 
rest of the stellar surface. We therefore consider the specific 
intensity of the spot covered surface to be described by limb 
darkening coefficients di, d^, dz and d^: 



C2.2 Case MS 

For case M3, it may be shown fsee lKippindl20l"ir ) that: 
Ftotai = / 2rl„(r)dr, 

Jr=0 

= i-Y^-. (en) 

^ n + 4 v ' 

n— 1 

For the spot, we must compute the flux obscured by 
its presence. This can be done by exploiting of the circular 
symmetry of the limb darkening effect and integrating the 
flux over an annulus defined to have an inner radius equal 
to the left-most point of the spot and outer radius equal to 
the right-most point of the spot. 

/•sin(/3 + ct) 

^obscured, annulus / 2r7* (T) dr, 

J r=sin(/3 — a) 

- cos 11 * 1 (13 + a)j . (C12) 

Note that the solution above has a significan tly more 
compa ct form that than acquired for an exomoon in lKippingl 
|201lf ). We may correct for the fact this is the flux over the 
entire annulus by applying: 
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Table CI. List of cases identified bv lMandel fc Agoll (|2002h . We use the same case classification in this work, but altering the notation. 



Case 


Analogous Condition for a Planet 


Condition for a Spot 


Range 


Ml 


1 + p < Sp t < oo 


Case IV 


it/2 + a < P < 7T 


M2 


1 - p < S P , <l+p 


Cases II & III 


tt/2 - a < p < tt/2 + a 


M3 


p < 5p» < 1 - v 


Case I 


a < p < it/2 - a 


M9 


<S P , <p 


Case I 


< P < a 



p M3 

* obscured 



ttA 



A p M3 

r obscured , annulu 



A 



r old \ o/zDi M obscured, annulus i 

TT [cos 2 (p — a) — cos J (p + QjJ 

(C13) 



where A and F, 



obscured, annul 



us have been previously de- 



fined. 

Finally, we need to compute i^pot, the flux from the 
spot itself. As will be the situation for all cases, the deriva- 
tion for F^pot is precisely the same as that as was done for 
^obscured except that {ci, C2, C3, C4} — ► {di , c?2 , d,3 , di\ and 
we multiply the expression by /spot to account for the tem- 
perature difference: 



Correcting for the expanded annulus area, we find: 



F, 



A 



obscured 



TV A annul us 

A 



7rcos 2 (/3 — a) 



p M2 

-F obscured, annulus i 



P M2 

* obscure 



(C18) 



C2.5 Case Ml 

The final case, and the simplest, is when the spot is out-of- 
view, analogous to the out-of-transit planet. Here, we have: 



F M1 - n 

^-obscured 



(C19) 



rpMx _ r y rpMx 

-^spot — ./spot UUl x* obscured j 

C XI 



(C14) 



where the x label emphasises that it is valid for all cases. 



C2.3 Case M9 

M9 considers the case when the spot overlaps with the centre 
of the stellar disc. Here, we must adjust the integration limits 
of the annulus flux since r > 0: 



obscured , annulus 



sin(/3 + c<) 



r=0 
4 



E \4 + n) 



2rl, (r) dr, 



4c n \ r 

1 — cos 2 



(C15) 



Correcting for the expanded annulus area, we find: 



& pA/9 

r obscured, annulus i 



^"^annulus 

A 



tt[1 - cos 2 (/3 + a)] 



F 



obscured, annulus ■ 



(C16) 



C2.4 Case M2 

M2 considers the case when the spot now hits the stellar 
limb. Again, we must adjust the integration limits of the 
annulus flux since r < 1: 



C2.6 Final Expressions 

For a single rotating starspot, satisfying the assumptions 
made in this work, we find that the flux from a star with a 
starspot may be written as: 



E 



n+4 n + 4 

4(Cn - dn/spot) C- 2 -C+ 2 



n + 4 



where 



and 



C4 



if < p < a, 
cos(P - a) if a < /3 < § + a, 
,0 iff+a</3<7r, 



cos(P + a) if < P < f - a, 
if f - a < P < TT, 



(C20) 



(C21) 



(C22) 



In the above form, the expressions span two/three do- 
mains. A single-domain function can be expressed using 
Heaviside Theta functions, H(x): 



F M2 

r obscured , annulus 



sin(/3 — ol) 



2rh(r) dr, 



r=l 
4 



K*^-4 (on) 



C- = cos(/3 - a)H(p - a)H(- - (/3 - a)) + H(-(/3 - a)), 
C+ = cos(p + aW/3 + a)H(~ - (/? + a)) + H(-(/3 + a)). 



(C23) 



Or more generally: 
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or explicitly 



£(x) = cos a;H(a;)H(— — x) 



+ H(- 



(C24) 



where C- = C(P - a) and ( + = £(/3 + a). 

Equation IC20l may be shown to return Equation IC5I if 
one sets {ci, C2, C3, ca\ t = {di, t?2, ^3, di} T = {0, 0, 0, 0} T , as 
expected. For Equation IC51 we showed how it was trivial to 
generalise the expression to Ns spots, provided one assumes 
the spots do not overlap. The same extension may be used 
here to yield: 



fHfl-i-tGSi)-^ 



4:(C n dj n j 'spot) 



n+4 n+4 

C 2 — f 2 



n + 4 



fcjC— ,k 



(C25) 

where the expressions for C+/- are trivially generalized 
to by amending a — » a*, and /3 — > Note that in the 

above expression we have added a Kronecker Delta function. 
This is because for f3 > n /2 + a, the fraction containing the 
£ terms goes to 0/0 i.e. undefined. Adding the Kronecker 
delta instead causes this to be equal to 0/1 = in this 
special circumstance and thus adds numerical stability to 
the function. 



C3 Expressing /3 with Physical Parameters 

C3.1 Accounting for the Star's Geometry 

So far, we have derived an expression for the flux from a 
limb-darkened star covered in multiple spots of sizes a and 
instantaneous positions f3, as given in Equation IC25I It was 
shown earlier how /3 could be related to a specific choice of 
apparent longitude, A, and apparent latitude, $, via Equa- 
tion [MB 

P(k, l>) = cos -1 [cos Acosl>]. 

As stressed throughout, A and <3> are the apparent lon- 
gitude and latitude of a starspot. The vector describing the 
Cartesian coordinates of the spot's centre is R CC ntrc and so 
far we have only defined this as a function of A and $ i.e. 
we know R cen tre(A, However, here we show how the vec- 
tor can also be expressed as a function of the true longi- 
tude and latitude (i.e. accounting for the star's geometry), 
Rcentre(A, $). This is crucial since the flux from the star is 
described by the parameters a and /3 only and ultimately 
one wishes to describe the flux as a function of the physical 
parameters and not auxiliary angles. 

Consider a frame in which the geometry of the star is 
such that the rotation axis has a normal vector given by j i.e. 
along the K-axis. In this frame, the apparent longitude and 
latitude are in fact equal to the true longitude and latitude, 
by virtue of definition. In this frame, which does not account 
for stellar geometry, we describe the position vector of the 
spot's centre with vector R centre . Due to the argument made 
above, we have: 



Rcentre (A, $) = Rcentre (A = A, 3> = $), 



(C26) 





rv 

■''centre 




sin A cos $ 


■^centre 


/ 

2/ccntrc 




sin $ 




/ 

_ ^centre _ 




cos A cos $ 



(C27) 



In order to calculate Rcentre (A, "!>), we must transform 
the frame to account for the geometry of the star. In other 



words we seek to transform R' 



Rc 



Euler's rotation theorem states that any large series 
of three-dimensional rotations can be written as a series 
of just three rotations only. Two conventions exist for how 
these three "Euler rotations" may be performed. The first 
is known as "proper Euler angles'. According to the intrin- 
sic/extrinsic rotation equivalences, proper Euler angles are 
equivalent to three combined rotations repeating exactly one 
axis e.g. X-Z-X. The second convention is called the "Tait- 
Bryan angles" (also known as the "Cardan angles") and 
these are equivalent to three composed rotations in different 
axes e.g. X-Z-Y. 

We abstain from choosing a convention for the mo- 
ment and proceed to consider a sequential choice of rota- 
tions which minimises the degeneracy between the various 
angles. We note that by choosing the first axis to be Y, we 
can eliminate a redundant angle since we defined an initial 
configuration with the stellar rotation axis aligned to the 
K-axis (i.e. an initial Y rotation is equivalent to intrinsic 
stellar rotation). For the sake of completeness, we refer to 
this first rotation as a clockwise rotation about the Y-axis 
by an angle uj*. 

For the next rotation, it is desirable to include stellar 
inclination at this point. A clockwise rotation about X by 
an angle (n/2 — J*) would correspond to the traditional def- 
inition of the stellar inclination angle. We have now selected 
the first two rotations, leaving just to the third. If we follow 
the proper Euler angles, we will be forced to us a Y rota- 
tion. In contrast, the Tait-Bryan convention would require 
a Z rotation. Since the observer is located down the Z-axis, 
a rotation about this axis cannot change the observed disk- 
integrated flux. Thus, a rotation about this axis would be 
redundant. For this reason, we use the Tait-Bryan conven- 
tion and define our Euler rotation scheme as Y-X-Z leading 
to two redundant angles and only one angle of physical in- 
terest, 7* (for completeness we dub the Z rotation angle as 
ifi*). We therefore define the position of the starspot centre, 
after applying the Tait-Bryan rotations, as: 



R c 



c 



(C28) 



where the first rotation is a clockwise rotation about 
the Y-axis by an angle w„: 



cos ui* sin ui t 

1 
— sin cos uu* 



(C29) 



The second rotation is a clockwise rotation about the 
X axis by an angle (it/ 2 — 1*). 



M/, = 



1 

sin /* — cos /, 
cos /* sin 7* 



(C30) 
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Finally, the third Euler rotation is about the Z-axis in 
a clockwise sense by an angle 



M 



cos V>* — sin ip t 
sin ip* cos V>* 
1 



(C31) 



Recall from Equation IA21I that the Z-component of 
R-ccntro directly yields f3, via 



I -^-centre I 



(C32) 



Applying all three rotations and extracting the Z- 
component allows us to write j3 as a function of the true 
longitude and latitude: 

ft = cos~ 1 [cos(A + to*) cos <E> sin/* + cos/» sin<!>]. (C33) 

As discussed earlier, and manifestly evident from the 
above expression, the angle u, is fully degenerate with A 
and thus may be neglected, giving us: 

/3(A, $,/,) = cos -1 [cos A cos $ sin/, + cos /* sin $] . (C34) 

It may be easily seen that this is p recisely the same 
expression as Equation 8 of lDorrenl (ll987l V 



C3.2 Accounting for the Star's Rotation 

Stellar rotation causes the a spot's instantaneous longitude 
to vary as a function of time. We denote the rotation rate by 
f2» = 2ir/P t , where P, is the rotational period of the star. 
Although it may be possible to envisage spots which migrate 
in latitude as well as longitude, we here only consider the 
simple case of <j> = d$/dt = and A = dA/dt = fi,. We 
may then decsribe the spot's instantaneous longitude and 
latitude as a function of time using: 



A(t) = A(t = i re f) + A(t - tref,h) = A(t = tref) + 
2n(t - trrf) 



A ref + 



P, 



2Tv(t — tref) 
P 1 

(C35) 



$(t) = $(t = tref) + $(t - tref.fe) = $(t = tref), 

= $ref. (C36) 



C4 Differential Rotation 

In general, unique P„^ terms are included to allow for 
differential rotation. However, differential rotation is well- 
described by the following function: 



P*,k = 



1 — K2 sin $ r 



K4Sin $rcf,fc 



(C37) 



The sin 4 ^ r d,k is usually only used for Solar rotation 
analysis, but we include it here for cases where a user re- 
quires a more sophisticated differential rotation profile. 



C5 Starspot Evolution 

C5. 1 A Linear Model 

There is some debate in the literature regarding 
how to model the evolu tion of a single starspot. 
iRiidiger fc Kitchatinovl (2000) show that the theoretical de- 
cay rate from 2-D modelling of a sunspot is close to lin- 
ear fo r the spot a rea (i.e. a square-root rate for a). Sim- 
ilarly, IStixl (|2002h argue that if the decay of an isolated 
sunspot is set by the amount of the azimuthal electric cur- 
rent with in the spot, a linear decay in a r ea wo uld result. 
However, iPetrovav fc Van Driel-Gesztelvil l|l997f l use a sta- 
tistical analysis of sunspot data to show that an "idealised" 
sunspo t exhibits parabolic are a decay. Additionally, observa- 
tions of lMartmez et af] (|l993h find both linear and parabolic 
decays. 

With our model, the user is free to use any descrip- 
tion they desire, but the code provided considers a simple 
linear growth/decay in a model and the partial derivatives 
computed are only valid for said model. 

A linear model has the advantage of being intuitively 
simple to handle, making the selection of appropriate bound- 
ary conditions and starting points easier for regression prob- 
lems. It is also very quick to computationally evaluate and 
encapsulates the key physics involved. Our linear model pro- 
duces a linear growth, flat-top and then linear-decay i.e. a 
trapezoidal profile for the spot's evolution. The profile is 
allowed to asymmetric to reproduce realistic evolution. 

We model starspot growth/decay via the a parameter 
only i.e. we consider the flux contrast to be constant. The 
linear-model has the simple form: 



Qfc(ti 

f^max.fc 



= I t T 1 [AtiH(At 1 ) - At 2 H(Af 2 )] 
- £: fc " 1 [At 3 H(At 3 ) - At 4 H(At 4 )]. 



and using 



Atl — ti tmax,fc 

At 2 = U - 

tmax,fc 
A^4 = ti tmax,fc 



~~t +lk ' 



,k + —, 

_ Lk 
,h 2 , 

Lk 
2 



(C38) 

(C39) 
(C40) 
(C41) 
(C42) 



where H(s) is the Heaviside step function, a max ,fc is the 
maximum spot-size, Lk is the full- width-full-maximum "life- 
time" of the spot and Tk & £ k are the ingress & egress du- 
rations of the spot profile. 



C6 Normalisation and Blended Light 

As discussed earlier, we choose to normalise the flux of spot- 
ted star by the flux which the same star would cause if no 
spots were present i.e. 

_ F(a,/3) 

Fmod - F( a = o,py (C43) 

At the time of writing, the most precise and sizable 
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source of photometric time series for main sequence stars 
comes from Kepler Mission. With this in mind, we choose to 
include a blending factor at this stage to account for overlap- 
ping PSFs, background flux or flux from an associated mem- 
ber in the system. This is a fairly common occurrence for 
Kepler data and many other photometric surveys due to the 
crowded fields observed. Let us consider then that the total 
flux observed is changed via F(a,f3) — > (F(a,0) + Fbiend). 
Our normalisation factor must now also be modified if we 
require that F mo d = 1 for an unspotted star. An appropriate 
choice is to use {F(a — 0,/3) + Fbiend): 



F(a,0) + Fbiend 
F(a =0,/3) + Fbiend ' 



(C44) 



Using Fbiend is cumbersome and a more common ap- 
proach is to define a blendin g facto r, relative to the tar- 
get's flux. iKipping fc Tinettil (|2010h advocate using B = 
(Source + fbi en d ) /Source which we follow here. This yields: 



F(a,0) 



B-l 



B 



BF(a =0,0) B 
F(a = 0,0) + fbiend 
F(a = O,0) 



(C45) 
(C46) 



CT Allowing for Multiple Time Series 

In many practical cases, we must fit multiple epochs of data 
which have different systematics. The most common sys- 
tematic to be treated is a baseline parameter, U. A common 
application of this process is using a unique normalisation 
factor for each Kepler quarter since spacecraft rolls affect 
the total flux within a defined aperture. For M data sets 
(e.g. M quarters of data from Kepler), each set requires a 
unique U m parameter. Using a box-car function (II), which 
is a composite of two Heaviside Theta functions, one can 
reproduce the desired behaviour: 



M / 

u m n m 

m=l \ 



F{oc,0) + B„ 



B m F(a = O,0) B„ 



, (C47) 



T ejx d,m) = H(t — T £ 



start, m ) H (t Tcnd.m ) ? 

(048) 

where it is understood that T sta rt,m+i ^ Tend,™ > 
Tstart.m- Note that we have also assumed that each data set 
has a unique blending factor. For Kepler data, it is typical 
for each quarter to have a unique B factor from spacecraft 
motion altering the PSF overlaps. 

Consider we have two data sets separated by N rota- 
tion periods where N 1. Further assume that the time 
span of data sets 1 and 2 are shorter than the spot life- 
time of all spots i.e. (T cnd ,m — Tstart.m) < Lk for all k 
and m. In this case, may one wish to treat the spots in 
data set 1 as independent of data set 2. This can be im- 
plemented by making use of the starspot evolution equa- 
tions. Specifically, one wishes to impose box-car spots (un- 
changing during each data set) with cut-offs in-between 
the two data sets. So the k th starspot of the ro' data 
set will have a m ,fc(t; a m ax,m,k, t m ax,m,fe, im,t, s mjfe ) take the 
form s m>fe — > co and i max , mifc = (T cnd , m — T Btalt ,m)/2 and 

Lm,k — (Tcnd,m — Tstart ,m) ■ 



APPENDIX D: PARTIAL DERIVATIVES 
Dl Motivation 

One of the major benefits of writing our expression as a 
single-domain function is that one can consider writing down 
a set of a single-domain partial derivatives. Partial deriva- 
tives are highly useful in optimisation problems where fre- 
quently the Jacobian matrix is computed to expedite a re- 
gression problem. 

D2 Partial Derivatives of the Likelihood Function 

The commonly used Gaussian noise likelihood has the form: 

c(&) = ft -jL= e* P [ - (F ° bs - - F r- (e))2 i . (di) 

i=1 V 27rcr i L M i 1 

Taking the partial derivative of the log likelihood with 
respect to parameter Qj yields: 



dlog£ 



i = l 



d(ri/ai) 



dF m 



(99, 



(D2) 



where ri = (f bs,i — fmod.O- Also note that in the above, 
and what follows throughout, that any partial derivatives 
taken with respect to Qj implicitly means that all other 
parameters are held constant except Qj . In other words, for 
a set of parameters O; where I = 1,2, ...L — 1,L, we use 
the notation that the partial derivative of some quantity X 
follows 



OX 
<99, 



dX 
<99, 



(D3) 



The outstanding problem is to derive dF mo d,i/dQj, 
which we deal with in the next subsection. We point out that 
any reasonable likelihood function, even if non-Gaussian, 
will still require dF mo d,i/dQj. For this reason, in the pro- 
vided code macula, we do not provide the partial deriva- 
tives of a Gaussian likelihood function directly but instead 
provide the partial derivatives of the model flux instead, 
9fmod,i /dQj . In this way, the results from macula are more 
general and hopefully of greater use to typical observers. 

D3 Partial Derivatives of the Model Flux 

Recall our final expression for the model flux, evaluated for 
the i th data point: 

_ A ( F( a ,(3) Bm-l \ 

- 2^ y BmFi{a = o, & + Bm )' 



which may be written as 

M 

fmod.i ^ ^ F n 



(D4) 



(D5) 
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So one may easily see that 



dF m 



E 



<9F mod 



(D6) 



The Fmod, m ,i function now requires partial derivatives. 
We adopt the assumption that II m ,; is not a function of any 
of the Qj parameters. This is perfectly reasonable as the 
function is only a function of T s t ar t,m and T cn d,m, which the 
user would define rather than fit for. Using this assumption, 
dH m ,i/dOj — for all i,j,m. Using the replacement (purely 
to save space) that F 4 = Fi(a,f3) and Fo,, = F»(a = 0,/3), 
one may now show: 



dQ 



m = U B ™ + ¥ »AB m - 1) 

+ U m ( B m ¥ i — B m ¥ 

+ Fo,(Fo,-F ^)). 



dU m 



d¥ ,j 



(D7) 



The above expression shows that the partial derivatives 
of the model flux can be expressed as a function of four other 
partial derivatives (which in turn may be broken down into 
other partial derivatives). 

Fo,» and in particular F» arc functionally dependent 
upon many Qj parameters but U m and B m do not. Rather, 
they represent a fitted parameter and have no other depen- 
dencies. We therefore have: 



OUrr, 

<96, 



and 



8Q 



if Qj U m , 

1 if Qj = Um, 



if Qj^B m , 
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j | 1 if Qj = Br, 



(D9) 



With these expressions the only remaining partial 
derivatives to find are those of F; and Fo,;. In fact, since 
Fo,« is defined as simply a special case version of Fj then we 
only require solving the partial derivatives of Fj or equiva- 
lently F(a,(3). 



D4 Partial Derivatives of the Flux w.r.t. Limb 
Darkening 

The Fi(a,f3) function is fully expressed as: 



4 N S A 



4(c n — dn/spot,fc) 



C— , h , i ^+ , k , i 



" + 4 C-,fc,i-C+, fc ,i+*C + ,M.C-,M 



For an unspotted star, A kti = for all k & i and so one 
may write: 



F i (a,f3) = F i {a = 0,f3)-Qi 

«(«• .ffl-'-t^) 

n— 

N S 

Qi = ^2 qk,i 



(D10) 
(Dll) 



(D12) 



Qk,i — 



E 



4(c n G^n/spot,fc) 



" + 4 n + 4 

A 2 _ A 2 
"»-,*,» ^ + ,k,i 



(D13) 



This allows us to write that 
0Fi(a,/3) _ 9Fi(a = 0,/3) 



96, 

9^(a^0,/3) 
96, 



ae 3 

iV s 



E U<Jk,, 



(D14) 



It is easy to show that 



dF i (a = 0,/3) 
6, 



-5 if6j=ci, 

-5 if 63=02, 

-f if 6^=03, 

-i if 63=04, 

otherwise. 



(D15) 



D5 Partial Derivatives of q k ,i 

The outstanding problem is now to find the partial deriva- 
tives of qk,i with respect to 6j. q k ,i is defined as: 



A k ,i I ^""^ 4:{Cn ^n/spot,fc) 
TV \ 



n+4 n+4 

f 2 — c 2 

>—,k,i *+,k,i 



n + 4 



C+,fc,i,C-,fc, 



We therefore proceed to derive the full four-coefficient 
partial derivatives, which we start by re-writing: 



qk,i — > Wn.k.i 

IT 

n=0 

^{Cn d n f 'spot, k) 



(D16) 



, k , i ,k,i 



" + 4 C,fc,i-C^, M + *C+.M.C-,* 



(D17) 

For the complex function _4(a,/3), the only derivatives 
of interest are with respect to a and /3 since A is functionally 
dependent on these terms alone. It may easily shown and 
numerically verified that: 



dR[A] \ 
da ) 



d/3 



OA 

da 

OA 

0/3 



(D18) 



(D19) 



Since all other partial derivatives of A can be expressed 
using the chain rule as a combination of the above two forms, 
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then partial derivatives of A can be derived for all Qj using 
this simple trick. This allows us to write: 



J n=0 J J n=0 

(D20) 

With the above, one can see the outstanding problem 
is to find partial derivatives of Ak > & t»n k i with respect to 



D6 Partial Derivatives of Ak.i 

Ak,i is a function of a k ,i and f3 k ,i only. Whilst these two 
terms will be functions of other parameters themselves, they 
offer the obvious starting point for a derivation of Ak,i's 
partial derivatives. The partial derivatives with respect to 
a k ,i and fi k ,i are easily shown to be given by: 



dA k ,i 
da k ,i 



dAk,i 
d/3 k}i 



= — sino^i sin (3k,i£ k ,i + 2cosa k ,i cos Pk,i=-k,i, 

(D21) 

= 0.5cosa fei icos/3 fci ie fci i - sina^j sin/3 fc>i H fcii . 

(D22) 



where we use 



£fc,i = 



esc 2 p k ,i(cos 2a k ,i + cos 2/3k,i) 



(D23) 



One may now employ the chain rule to write: 
dAk, i 



+ 



dAk^ 




( da k ,i 


da k ,i j 




\ d®j 


dA k ,i\ 




( dp k ,i 


dp Ki J 


a k,i 





(D24) 



where we temporarily re-include the implicit notation 
to make the expression less ambiguous. Partial derivatives 
of a k ,i &i fi k ,i with respect to Qj will be provided later. 



D7 Partial Derivatives of w n , k ,i 

w n ,k,i is fully expressed as: 



Wn,k, 



T 



n,k,i 



4(C n d n /spot.fc ) 

71 + 4 

C 2 — C 2 , 

,k,i = + i 



(D25) 
(D26) 



We first turn our attention to taking the partial deriva- 
tives of Y n ,k,i with respect to Qj. We note that that: 



<><\ 



da k ,i 



which via the chain rule imply 

dS C+,k, t X-, k ,i 



o. 



= 



(D27) 
(D28) 



<99, 



V{i,j,k} 



(D29) 



With this simplification, we find: 



96, ^+' fc ' < 90,- 



2T n t j 



The partial derivatives of C-,fc,t are given by: 



(D30) 



dotk,i 



8[—{Pk,i — a k ,i)] 



Sk,i 



H(f - (13k,, -a k ,i)M/3k,i ~ak,i) 

csc(/3 fc ,i - a fc]i ) 
2% - 2(/3 fc ,j - Qfc,i)]H(/3 fc ,j - Q fc ,i) 

sec( / 9 fe , l - afc.i) 
5[/f3fc,i — a feii ]H(§ — (/3fc,i — Qfe.i)) 



sec(/3 M - a fc)i ) 



' «fc,« 
and of C+,fc, 



9C- 



<9a fc ; 



, (D31) 
(D32) 



dC+,k,i 

da k i 



Pk,< 



2S[tt - 2{p k ,i + Qfc.i)] 
sec(^ fc ,i + afc.i) 

H(/3 M +a M )H(f - (^.i + a fc ,i)) 



csc(/3 fcii + a fc]i ) 



^C+,k,» 



^C+,fc,i 

<9a fejl 



(D33) 
(D34) 



In practice, the 5(x) functions always yield zero unless 
x = 0. Since they are a function of a continuous variable, 
namely time, the probability that the time will precisely 
yield a non-zero S function is infinitesimal. For this purpose, 
they are simply ignored in the macula code. The latter rela- 
tions lead to a simplification of the chain rule expansion of 
the partial derivatives with respect to Qj : 



An Analytic Model for Rotational Modulations 27 



d(-,k,i 




de j 


y 9a k4 




( d(+,k,i 




i da k ,i 



13k,, L 



da k ,i _ dfik,% 

ae 7 - 



ae, 



da k ,i d/hj 

ae, ae, 

Finally, the partial derivatives of w n ,k,i are 



(D35) 
(D36) 



e, 



n + 4 



1 n,k,i 777 ~r (^Cn Ctn/spot.fc J " 



96, 



n J- n.kA 777 ,/ spot ,fc J- n,fc,z 7777 

O0j aQj 



ae, 



(D37) 



Since the partial derivatives of T n ,k,i have been dealt 
with above, this leaves us to comment on the partial deriva- 
tives of Cn, d n and / spo t.fc- These represent fitted parameters 
(or perhaps fixed) and thus one may trivially evaluate their 
derivatives to be 



dc„ 



if / c n , 

1 if 6, = Cn, 



(D38) 



d/3k,i _ sin $ rc f ,fe sin I, — cos A k ,i cos <£> rcf , fc cos 7* 
sin/3 M 



dh 

dp k ,i 
OPbq 

dK2 

a«4 

a/3 fc , t 
a<E> rcf ,fe 



(D42) 

27r(t» -t r cf,fc) cos<j> rc f !fc sinA M sin I, (t^aq) 

ft:Q-P.,fc sin/Jt,,; 
27r(ti - t rcf , fc ) sin 2 $ rcf , fc cos$ rc f,fc sin A M sin J, 
-Peq sin/3 fc ,i 

(D44) 

27r(tj — t rc f,fc) sin 4 $ rc f,fc cos <£ re f,fc sin A*,,, sin J, 
Peq sin/3 fe>i 

(D45) 

= esc /3 k ,i sin /» sin <I> rc f , fc ( cos A fc>i 



sinA fe]i [2K 2 cos $rcf,t + K4Sin (24> rcfifc ; 

^EQ 

- csc/3 fcii cos h cos$ re f,fc, (D46) 



d/3 k<i _ sin/, cos <i> rcf , fc sin A fc|i 



. (D47) 

aA rcfife sm/3 fc ,i 

above, the remainder of the partial 



Aside from the 
derivatives satisfy: 



= if 9 3 + I,,P E Q,«;2,K4,$rcf, fc ,A rcfifc . (D48) 



dd IL _\o ife./dn, 

dQj ~\i if 9j ■ = d„, 



D9 Partial Derivatives of a*,,, 

(D39) The k th starspot evolves via Equation [5] The partial deriva- 
tives are found to be: 



a_/spot ,k 

ae, 



if Oj 7^ /spot,*:, 

1 if Bj = /spot,*- 



(D40) 



D8 Partial Derivatives of 

The only partial derivatives now missing are those of ctk,i 
and /3 k ^ with respect to the fitted parameters, 0j. /3 kti is 
defined as: 



Pk,i — cos 1 |^cos/* sin$ fc ,i + sin J» cos$ fc ,i cosA feji j , 

A _ a 2ir[ti — tref.fe) 

Afe,i — A Ie f t k H 75 1 

Accounting for differential rotation, the longitude evo- 
lution is described by: 



* . , 2iv(ti — t Te !,k) ,-. . 2 ■ 4 « \ 

Afci = A re f ,fc "I (1 - K2 SU1 $rcf ,fc - K 4 Sin <P re f ,fc ) 



EQ 



(D41) 



Now the partial derivatives yield: 



Q-k,- 



da k ,j 

^max , k Oim&x, k 

da k ,i 



dt n 



2s»!*(H(A*i) - H(At 2 )) 
+ 2w ! * (H(Ats) _ H(At4))j 

^ = r w^ (H (At 1 )-H(Afa)) 



(D49) 



(D50) 



+ 



da k 



dl k 

da k ,j 
d£k 



H(Ai 3 ) - H(Att)), 

H(Ati) 



2£ k 

a ma x,fc(Ati + At 2 



2Z 2 
fe (At 3 + Ai 4 ) 



(D51) 

H(Afe)), 
(D52) 



2f 2 



;H(Ai 3 ) - H(At*)). 



(D53) 
(D54) 



Aside from the above, the remainder of the partial 
derivatives satisfy: 



da k 

ae. 



if e, + 



(D55) 



Finally, it is necessary to define a reference longitude, 
A re f, k- A convenient choice is to define it as the longitude 
at the instant t — t maXj fc, which is the default assumption of 
macula. 
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APPENDIX E: MANIPULATING PARTIAL 
DERIVATIVES 

El Alternative Limb Darkening Laws 

In this model, we have adopt e d the four-coefficient limb 
darkening proposed by I Claret! l)2000l ). Several, but by no 
means all, alternative limb darkening laws can be adopted 
by re-parametrising the four-coefficient model presented in 
this work. In this subsection, w e discuss four exa mples: i) 
the q uadratic law [|Kopa]|ll950h ii ) the linear law ijRusselll 
Il912h iii) the thr ee-coefficient law (ISing et alj|200gh iv) the 
square-root law (jDfaz-Cordoves fc Gimeneall992i ). In each 
case, we show how one may use the results from macula 
to obtain the partial derivatives for regression purposes. In 
what follows, we assume that the spot and the star have 
distinct limb darkening coefficients. 



< m < 2, 
< u 1+2 < 1. (E4) 

The four-coefficient model can be set to these parame- 
ters using: 

ci = 0, 

C2 = 2ui+2 — Ml, 
C3 = 0, 

C4 = Hi — Ul+2- (E5) 

In the previous section, we have derived dF mo d/ dc n for 
n = 1, 2, 3, 4. We now require dF mo d/du\ and <9i ? mo d/<9iii+2. 
Firstly, one may show: 



El.l Quadratic Law 

The quadratic law, first proposed by iKopall (|l950h . is per- 
haps the most commonly adopted model in the exoplanet 
literature. Recall from Equation [10] that the four-coefficient 
limb darkening law is described by 



/»(r) = l-^> n (l-^ 2 ). 

71=1 

In contrast, the quadratic law is described by 



/„(r) = l-Ui(l-M)-W2(l-M) 2 - (El) 

By comparing the coefficients relative the four- 
coefficient model, one may show that the quadratic law may 
be reproduced by setting: 



ci = 0, 

C2 = Ml + 2ll2, 
C3 = 0, 
C4 = — U2. 



(E2) 

The quadratic model is popular in the exoplanet com- 
munity when one wishes to fit for the limb darkening param- 
eters. The reason for this is two-fold. Firstly, photometric 
data are rarely precise enough to regress a unique solution 
for all four coefficients of the non-linear limb darkening law 
and so using the quadratic model reduces the number of 
free parameters by two yet preserves curvature in the in- 
tensity profile of the star. Secondly, the two quadratic coef- 
ficients, u\ and it2, have well-described priors by imposing 
that the intensity profile is monotonic and everywhere posi- 
tive. ICarter et all (2009) show that these conditions impose 



ui = c 2 + 2c 4 , 

Ul+2 = C 2 + C 4 . 

It is therefore trivial to write: 



(E6) 



dui 

9-Fmod 



-\- Z- 



dui 



dC2 
dF mod 
dC2 



+ 



<9C4 
dF mod 
dC4 



(E7) 



For the starspot's limb darkening, the same argument 
may be made to show: 



dF a 



dvi 

9-Fmod 



9F mo( \ ~ 0F mo d 

■ -p Z - 



dd 2 
dF mo d 



ddi 
dF mo d 



dv 1+2 
where we define di 



dd2 ddi 
d-j, — and 



(E8) 



Vi = di + 2d 4 , 
v 1+2 — d 2 + di. 



(E9) 



El. 2 Linear Law 



The lin e ar lim b darkening law, which can be traced back to 
iRusselll l)l912h . is expressed as: 



I,(r) = 1- Ui (l-M). 



(E10) 



It is therefore trivial to see that this is identical to the 
quadratic law where tti = ul and 112 = 0. Relative to the 
four coefficient model, we have {ci, C2, C3, C4} = {0, ul, 0, 0}. 
In such a model then, one may simply use: 



Ml > 0, 

0<mi+m 2 <1. (E3) 

iKipping et all (|2012h point out that a sensible upper- 
bound on mi may be imposed from insp ection of typ ical 
coefficient tables presented in works such as IClaretl (|2000l l. A 
typical choice is Mi < 2 for Sun-like stars. With this upper- 
bound one may re-define M1+2 = mi + M2 and regress the 
parameters {ui,mi+2} subject to the uniform priors: 



dF a 



dF^ 



dui 



dC2 



(Ell) 



As before, this can be easily applied to the starspot's 
limb darkening too via 



9F mo d 9F mo d 

du L dc 2 
where we define vl = (fe and d\ = d3 = d± = 0. 



(E12) 
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El. 3 Three-Coefficient Law 

The three-coefficient law, proposed by ISing et alj |2009l ). is 
described by: 



c 2 (l - n) - c 3 (l - m 3/2 ) - c 4 (l - m 2 ), (E13) 



which is precisely the same as the four-coefficient law 
in the limit c\ — > 0. For this reason, the partial derivatives 
are unchanged from before and one may ignore the partial 
derivative with respect to c\ & di- 



El. 4 Square-Root Law 

The three-coefficient law, proposed 

iDfaz-Cordoves fc Gimened Ijimk . is described by: 



h{r) 



ci(l - fi 



l/2x 



C 2 (l 



by 



(E14) 



which is again identical to the four-coefficient law in 
the limit of C3 — > and C4 — > 0. The same applies, of course, 
for the spot's limb darkening profile and thus the partial 
derivatives are trivially obtained. 

E2 Allowing for Global Parameters 

E2. 1 Principle 

Practically speaking, it is common to consider a subset of 
the parameters to be equal to some global term. For ex- 
ample, rather than regressing for Ns unique spot contrast 
fluxes, / S pot,fc, one may wish to enforce the condition that 
all spots have the same flux contrast (i.e. temperature). The 
advantage of implementing such a condition is that one re- 
duces the number of free parameters in the regression by 
{Ns ~ I)- 

In such a case, one requires the partial derivatives of 
the model flux with respect to this new global parameter, 
rather than the individual terms. Since the individual par- 
tial derivatives have already been calculated and are directly 
returned by the macula code, it is highly advantageous if we 
can phrase the partial derivatives of this new global param- 
eter as a function of the individual terms. In this subsection, 
we provide a method for accomplishing this. 

The model flux is a function of parameters i.e. 
i ? mod(0). For L model parameters, one may write out the 
differential as: 



And so by equivalence of Equations IE15KjEl6l one can 
see that: 



dG 

(=1 

And finally this yields: 

9-Fmod ,. 

= hm 

dG Si.ej Sy-jG 

E2.2 Common Examples 



dG, 



(E17) 



E 



8F a 



(Ei8) 



As we cited earlier, a common application of Equation IE18I 
is to Ns individual spot contrast values, f spo t,k to be equal 
to some global spot contrast term, g sp ot- The partial deriva- 
tive of the model flux with respect to this new global flux 
contrast new may be expressed, using Equation IE18I as: 



dF m 



lim 

/spot.k - ►ffspot 



E 



(E19) 



Another example is to enforce a global blending factor, 
C, rather than individual values, B m : 



<9Fmod 

dC 



lim 

B m ->C 



E 



dFrr 



8B n 



(E20) 



Finally, one may wish to set the spot's limb darkening 
parameters to be equal to the star's limb darkening param- 
eters i.e. c = d = b where b is the global limb darkening 
parameters in vector-form. 



— — = lim lim 



dc„ 



(E21) 



APPENDIX F: PARTIAL DERIVATIVES WITH 
RESPECT TO TIME 



Recall from Equation ID4I and Equation ID5I that the m 
component of the model flux is given by 



(=1 
dF mod 
90i 



dQi 



<56i 



dF„_ 



d<d 2 



-80 2 + 



8F a 



-SQ L . (E15) 



Now consider that a subset of the © parameter vector is 
set to be equal to some global parameter, G. Let this subset 
run from parameter 1 to L' i.e. Oi — O2 = ••■ = ©l = 
G where G is some global parameter. The differential now 
becomes: 



SFn 



dF„ 



dG 



lSG + E ^||f 5Qi - (Ei6) 



,. . t/ m n m ,jFi (B m — i)iy m n m ,i \ 

rmod,ra,i = I ^ — ^ 1 ^ • (r 1J 



BmFo 



Taking the partial derivative of the above with respect 
to time yields 



8F, 



mod,m,i 

~dti 



UrrJl m ,i 9¥i 

B m W ~dt~. 



+ 



u m w t Um(B m - 1) ^ an m ,i 



B„ 



dt, 



(F2) 



The partial derivatives of the box-car function, H m ,i, is 
simply two Dirac Delta functions and thus may be neglected 
in what follows i.e. 
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CP3) 

dU B m F dti' y ' 

Since Fj = Fo — Q% and Fo has no time dependency, 
then dWi/dti = —dQi/dU giving 



<9Fi _ _ dQk.i 



dt, 



(F4) 



The qk,i partial derivative may be expressed via 

Ah,, 



Qk,i 



^ W n ,k,i, 



dq k ,t _ M.[Ak.i] dw n>k ,i 1 \dA k ,i 



dU 



+ —i 



n— n— 



4 



(F5) 

For the partial derivatives of Ak,i, we can use the same 
chain rule trick as was used earlier: 



The partial derivatives of C-,k,i and C+,k,i with respect 
to ak.i have already been calculated in Equation ID32I and 
Equation ID34I respectively. The outstanding task is now to 
compute the partial derivatives of cikj and /3k,i with respect 
to time. It is easy to show that the ak,i partial derivative is 
given by: 



dak,-. 



Q^max, fc 

Ik 



(H(Ati) - H(Ata)) 
(H(At 3 ) - H(At 4 )), 



(Fll) 



and that of fik,i by: 

dfik.i _ 2m sin I, cos <j> r cf,fc sinA M 

dU 



P*, k yl — (sin /„ cos <f> TC f,k cos K k ,i + cos 7* sin $ rc f,fc 

(F12) 



This paper has been typeset from a Tp^X/ I^TgX file prepared 
by the author. 



dA k .j 
dU 



e, v! 



+ 



OAk.r 

dctk,i 
dA k ,i 



da k ,i 
dU 

dfik.i 
dU 



(F6) 



where the partial derivatives of Ak,i with respect to 
ctk,i and are given in Equations ID22I Let us leave aside 
the issue of the partial derivatives of aik,i and fik,i for the 
moment and focus on those of w n ,k,i'- 

_ 4(c n — rfn/spot.fc) ^ 

^n, k ,i , . J- n.fe.ij 

n + 4 

dw„ t k,i _ 4(c n — d n /spot,fc) d~£ n ,k,i Cw7"\ 
9ti ~~ n + 4 dti ' 

Partial derivatives of T n ,h,i with respect to 0; have al- 
ready been calculated earlier in Equation ID301 in terms of 
the derivatives of C-,fc,i an d C+,M- This result is easily mod- 
ified to be with respect to time: 



dT n , k 
~dt~ 



- ,k,i 



^ — , k , i ^j-/. ^ + , k , i 



2T r 



2 / 0^ 



dU 



dt, 



(F8) 



Those terms have also had their partial derivatives com- 
puted wth respect to ak,i and Pk,i, which lead to the chain 
rule relation: 



9C-,k,i 


( dC-,k,i 


dU 


I da k ,i 


9C+,k,i 


( 9C+,k,i 


dU 


\ da k ,i 



dce k ,i _ dp k ,i 
dti 

da k ,i 



dU 



+ 



dti 
dPk, 



dt, 



(F9) 
(F10) 



